#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <math.h>  
#include <map>
#include <fftw3.h>
#include <vector>
#include <iostream>
#include <iomanip>
#include <string.h>
#include <fstream>
#include <sstream>
#include <stdexcept>
#include <cstdlib>

#include <unistd.h> 
#include <stdarg.h> 
#include <errno.h> 
#include <fcntl.h> 
#include <time.h> 
#include <sys/time.h> 
#include <sys/mman.h> 
#include <sys/ioctl.h> 
#include <sys/stat.h> 
#include <sys/types.h> 
#include <sys/wait.h> 
#include <signal.h> 

// root includes
#include <TSystem.h>
#include <TROOT.h>
#include <TLatex.h>
#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <TBranch.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TProfile.h>
#include <TChain.h>
#include <TGraphErrors.h>
#include <TGraph.h>
#include <TGraphAsymmErrors.h>
#include <TEllipse.h>
#include <TApplication.h>
#include <TMultiGraph.h>
#include <TCanvas.h>
#include <TDirectory.h>
#include <TKey.h>
#include <TStyle.h>
#include <TPaveText.h>
#include <TText.h>
#include <TMath.h>
#include "TString.h"
#include "TVirtualFFT.h"
#include "TAxis.h"
#include "TPad.h"
#include <TPRegexp.h>
#include <TMinuit.h>
#include <TBasket.h>
#include <TLegend.h>
#include <TNtupleD.h>

#include "constants.h"
#include "inputtreevars.h"
#include "inputconfigvars.h"
#include "outputtreevars.h"

using namespace std;

void Put_Pointer_into_Vector(vector<double> *Vp, vector<double> & V);

void StripComments(string& );

bool IsSaturated(int fRunId, int fEvtId, int ich, vector<double> *ch, double fixvertgain);

bool SelTriggers(int fRunId, int fEvtId, double fFixVertGain[4], 
		 vector<double> *ch1,  vector<double> *t1, 
		 vector<double> *ch2,  vector<double> *t2, 
		 vector<double> *ch3,  vector<double> *t3, 
		 vector<double> *ch4,  vector<double> *t4,
		 double ne_min, double ne_max, double tb1, double tb2, double ts1, double ts2,
		 bool isdecember2012,
		 bool selch[5]
		 );

void TransformFFT(vector<double> *Vall, vector<double> *Tall, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq, vector<double> & F, vector<double> & P, vector<double> & VCut, vector<double> & T, vector<double> & PCut);

void TransformFFT2(vector<double> *Vall, vector<double> *Tall, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq, vector<double> & ReP, vector<double> & ImP);

void GetNe(vector<double> *sgn, vector<double> *time, double tped1, double tped2, double t1, double t2, double bunchlength, double &ne, double &ped, double &sgn_max, double &current);

double GetVrms(vector<double> *sgn, vector<double> *time, double deltaT, double t1, double t2);
double GetVrms1(vector<double> sgn, vector<double> time, double deltaT, double t1, double t2);
double GetVrms_bin(vector<double> *sgn, vector<double> *time, double deltaT, double t1, double t2); 

double PowerFromSpectrum(vector<double> *P, vector<double> *F, vector<double> MinFreq, vector<double> MaxFreq);
double PowerFromSpectrum1(vector<double> P, vector<double> F, vector<double> MinFreq, vector<double> MaxFreq);

void AlignSignals(double Tc1, double Tc2, vector<double> *fVA, vector<double> *fTA, vector<double> *fVB, vector<double> *fTB, double fNeA, double fNeB, double & fdt_alignment, double & fchi2);

double Aeff(double G, double nu);
double Calib_Dec2012(double nu);

double ExpectedVrms(double ne, double alpha, int antenna, double BunchLenght, int selfreq, vector<double> *time, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq);
double ExpectedVrms1(double ne, double alpha, int antenna, double BunchLenght, int selfreq, vector<double> time, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq);


TH1D *hh1_exp; TH1D *hh2_exp;
TH1D *h_AeffCal, *h_Aeff, *h_Cal;
double ExpectedVrms_MBRdEdX(double ne, double alpha, int antenna, int Targets, double BunchLenght, int selfreq, vector<double> *time, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq);


TF1 *fExpectedMBRVrms, *fExpectedMBRVrmsCut;
double GetCalibIntegral(int selfreq, vector<double> MinFreq, vector<double> MaxFreq);
double fMBRExpectedVRMS(double *x, double *par)
{
  
  double ne = x[0] * 1e9;

  double alpha          = par[0];
  int T                 = int(par[1]);
  double BeamLength     = par[2];
  double fCalibFactor   = par[3];
  double fRescaleFactor = par[4];

  // expected density flux
  double Ng = 1.20e7; // Gorham number of electrons.  
  double ratio = FrascatidEdXFoV[T]/SLACdEdXFoV[5];
   
  // Ghoram
  double I0 = 4.0*1e-16; // W/m2/Hz
  double Ls = 0.5;
  double Lf = 0.9;
  //double g = pow(Ls/Lf,2) * (1.3/3.0) * (10.0/BeamLength);
  double g = pow(Ls/Lf,2);
  double I = I0 * g * pow(ratio * ne/Ng , alpha);

  double P_exp = I * fCalibFactor * fRescaleFactor;
  return sqrt ( P_exp * 50 ); 

}


bool IsDecember2012 = false;

vector<double> MinFreq;
vector<double> MaxFreq;
vector<double> Fsgn, Psgn, VCutsgn, TCutsgn, PCutsgn;
vector<double> Fbkg, Pbkg, VCutbkg, TCutbkg, PCutbkg;


void InitOutputTree();
void BookOutputTree(TTree *fTreeOut);

bool DEBUG = false;


TChain fData("data");
// load the files and book the input tree
void BookDSTData();
/////////////////////////////////////////////////////

int main(int argc, char *argv[]) {
  //int main() {

  string dirname;
  vector <string> filename;
  filename.clear();


  cout << endl;
  ifstream fOpen(argv[1]);
  string buffer;
  TPRegexp isEmptyLine("^\\s*$"); // ie. only whitespace
  int icount = 0;
  while(!fOpen.eof()){
    icount++;

    getline(fOpen,buffer);
    StripComments(buffer);

    //Check if the buffer is empty
    if ( isEmptyLine.MatchB(TString(buffer.c_str())) ) 
      continue;

    //cout << buffer << endl;

    stringstream line(buffer);
    string descriptor;
    string value;

    //Get the name of the variable
    line >> descriptor; 

    // nota che   string dirname   deve essere messo prima del loop 
    // altrimenti in filename non me lo ritrovo !!!!! ?????
    // the name of the directory is always in the first line
    if (icount==1)
      dirname = descriptor;

    if (dirname.compare("~/AMY_DATA_Dec2012/")==0)
      dirname = "/Volumes/AMY/DSTFilesDec2012/";

    // get the DST file names 
    TString des1(descriptor);
    TString des2(des1(0,3));
    if (des2=="DST") {
      filename.push_back(dirname+descriptor);
    
      //Now check the date of the file
       TString tmp_date(des1(0,9));
       if(tmp_date == "DST201212")
	 IsDecember2012 = true;

    }


    TString des3(des1(0,2));
    if (des3=="#") 
      continue;

    // get the value of that variable!
    std::getline(line, value);
    TString buffer(value);
    // remove whitespace from beginning and end
    TPRegexp("^\\s+").Substitute(buffer, "");
    TPRegexp("\\s+$").Substitute(buffer, "");
    value = string(buffer.Data());

    //cout << descriptor << "\t" << value << endl;

    if (value.size()>0){ 
      if (descriptor == "BunchLenght[ns]") {BunchLength = atof(value.c_str()); cout << descriptor << "\t" << value << " --> " << BunchLength << endl;}
      if (descriptor == "Target") {Target = atoi(value.c_str()); cout << descriptor << "\t" << value << " --> " << Target << endl;}
      if (descriptor == "tmin_beam_ped") {tmin_beam_ped = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << tmin_beam_ped/ns << endl;}
      if (descriptor == "tmax_beam_ped") {tmax_beam_ped = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << tmax_beam_ped/ns << endl;}
      if (descriptor == "tmin_beam_sgn") {tmin_beam_sgn = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << tmin_beam_sgn/ns << endl;}
      if (descriptor == "tmax_beam_sgn") {tmax_beam_sgn = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << tmin_beam_sgn/ns << endl;}
      if (descriptor == "SelCh2") {SelCh[2] = atoi(value.c_str()); cout << descriptor << "\t" << value << " --> " << SelCh[2] << endl;}
      if (descriptor == "SelCh3") {SelCh[3] = atoi(value.c_str()); cout << descriptor << "\t" << value << " --> " << SelCh[3] << endl;}
      if (descriptor == "SelCh4") {SelCh[4] = atoi(value.c_str()); cout << descriptor << "\t" << value << " --> " << SelCh[4] << endl;}
      if (descriptor == "tmin_bkg") {tmin_bkg = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << tmin_bkg/ns << endl;}
      if (descriptor == "tmax_bkg") {tmax_bkg = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << tmax_bkg/ns << endl;}
      if (descriptor == "Ti0") {Ti0 = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << Ti0/ns << endl;}
      if (descriptor == "Tf0") {Tf0 = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << Tf0/ns << endl;}
      if (descriptor == "Ti1") {Ti1 = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << Ti1/ns << endl;}
      if (descriptor == "Tf1") {Tf1 = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << Tf1/ns << endl;}
      if (descriptor == "Ti2") {Ti2 = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << Ti2/ns << endl;}
      if (descriptor == "Tf2") {Tf2 = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << Tf2/ns << endl;}

      if (descriptor == "FreqCutType") {FreqCutType = atoi(value.c_str()); cout << descriptor << "\t" << value << " --> " << FreqCutType << endl;}
      if (descriptor == "FreqCutMin") {FreqCutMin = atof(value.c_str()); cout << descriptor << "\t" << value << " --> " << FreqCutMin << endl;}
      if (descriptor == "FreqCutMax") {FreqCutMax = atof(value.c_str()); cout << descriptor << "\t" << value << " --> " << FreqCutMax << endl;}

      if (descriptor == "Alignment") {Alignment = atoi(value.c_str()); cout << descriptor << "\t" << value << " --> " << Alignment << endl;}
      if (descriptor == "ich_ref") {ich_ref = atoi(value.c_str()); cout << descriptor << "\t" << value << " --> " << ich_ref << endl;}
      if (descriptor == "T_al_1") {T_al_1 = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << T_al_1/ns << endl;}
      if (descriptor == "T_al_2") {T_al_2 = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << T_al_2/ns << endl;}
      if (descriptor == "ReferenceSignal") {ReferenceSignal = atoi(value.c_str()); cout << descriptor << "\t" << value << " --> " << ReferenceSignal << endl;}
      if (descriptor == "Ti_all") {Ti_all = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << Ti_all/ns << endl;}
      if (descriptor == "Tf_all") {Tf_all = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << Tf_all/ns << endl;}
      if (descriptor == "Fmin_av") {Fmin_av = atof(value.c_str()); cout << descriptor << "\t" << value << " --> " << Fmin_av << endl;}
      if (descriptor == "Fmax_av") {Fmax_av = atof(value.c_str()); cout << descriptor << "\t" << value << " --> " << Fmax_av << endl;}
      if (descriptor == "tmin_plot") {tmin_plot = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << tmin_plot/ns << endl;}
      if (descriptor == "tmax_plot") {tmax_plot = atof(value.c_str()) * ns; cout << descriptor << "\t" << value << " --> " << tmax_plot/ns << endl;}
      if (descriptor == "Ne_norm") {Ne_norm = atof(value.c_str()); cout << descriptor << "\t" << value << " --> " << Ne_norm << endl;}
      if (descriptor == "Ne_min") {Ne_min = atof(value.c_str()); cout << descriptor << "\t" << value << " --> " << Ne_min << endl;}
      if (descriptor == "Ne_max") {Ne_max = atof(value.c_str()); cout << descriptor << "\t" << value << " --> " << Ne_max << endl;}

    }
  }
  cout << endl;
  fOpen.close();

  /**
  for(int ich=1; ich<Nch; ich++){
    if (SelCh[ich_ref] == false){
      printf("WARNING: SelCh set to true for the reference signal used to make the alignment \n");
      SelCh[ich_ref] = true;
    }
  }
  */

  string fin2 = argv[1];
  int index=fin2.length();
  fin2.replace(index-4, 5, ".root");
  TFile *fout = new TFile(fin2.c_str(),"recreate");

  // load the DST files 
  for (int j=0; j<filename.size(); j++){
    TString tfilename(filename.at(j));
    fData.Add(tfilename);
  }
  BookDSTData();

  TTree *fTreeOut =  new TTree("fTreeOut","fTreeOut");
  BookOutputTree(fTreeOut);

  //Load the calibration TGraph
  DoCalibration = true;
  TFile *fOpenCalibCables = new TFile("AttenuationTGraph.root","READ");
  TFile *fOpenCalibAntenna = new TFile("SatimoGainCalibration.root","READ");
  //TFile *fOpenCalibAntenna = new TFile("Horn.root","READ");
  TGraph *tgCable1 = (TGraph*)fOpenCalibCables->Get("File_R4");
  TGraph *tgCable2 = (TGraph*)fOpenCalibCables->Get("File_L5");
  TGraph *tgAMP = (TGraph*)fOpenCalibCables->Get("AMP41");
  TGraph *tgHorn = (TGraph*)fOpenCalibAntenna->Get("Horn120401A20");
  //TGraph *tgHorn = (TGraph*)fOpenCalibAntenna->Get("gain");
  tgCalib = new TGraph();
  tgAntenna = new TGraph();
  TGraph *tgAeff = new TGraph(); 
  for(int i=0;i<tgAMP->GetN();++i)
    {

      double freq, att_cable1;
      tgCable1->GetPoint(i,freq,att_cable1);
      
      double dummy_freq, att_cable2;
      tgCable2->GetPoint(i,dummy_freq,att_cable2);
      
      double att_amp;
      tgAMP->GetPoint(i,dummy_freq,att_amp);
      
      double antennaGain = 0;
      if(freq/1e9 < 17.99)
	antennaGain=tgHorn->Eval(freq/1e9);
      else
	antennaGain=tgHorn->Eval(18.);

      double tot_att = att_cable1 + att_cable2 + att_amp;
      
      //printf(" fHz = %.2f A4 = %.2f AMP=%.2f L5=%.2f AntennaG=%.2f\n",freq/1e9,att_cable1,att_amp,att_cable2,antennaGain);
      tgCalib->SetPoint(i,freq,tot_att); //dB
      tgAntenna->SetPoint(i,freq,antennaGain); //dB

      double area = Aeff(antennaGain,freq/GHz);
      tgAeff->SetPoint(i,freq,area);
    
    }
  
 

  fout->cd();



  ////////////////////////////////// time windows for power and Ne calculation //////////////////////////////////////////
  
  double tmin_sgn[Nch]  = {0., 0., Ti0, Ti0, Ti0}; double tmax_sgn[Nch]  = {0., 0., Tf0, Tf0, Tf0}; 
  double t1min_sgn[Nch] = {0., 0., Ti1, Ti1, Ti1}; double t1max_sgn[Nch] = {0., 0., Tf1, Tf1, Tf1}; 
  double t2min_sgn[Nch] = {0., 0., Ti2, Ti2, Ti2}; double t2max_sgn[Nch] = {0., 0., Tf2, Tf2, Tf2}; 

  //
  // oscilloscope traces in the selected time windows
  // il massimo intervallo temporale e il number trigger della traccia vengono definiti in alignment 
  // Se si effettua l'allineamento, il number trigger e' lo stesso di quello della traccia di riferimento usata per fare l'allineamento
  TGraph *gr_vt[Nch];
  TGraph *gr_vt_sel0[Nch]; TGraph *gr_vt_sel1[Nch]; TGraph *gr_vt_sel2[Nch];
  

  ////////////////////////////////////////////// Alignment in time /////////////////////////////////////////////////////

  //
  // The alignment is performed using Ch3
  // change by hand if you need to use another channel
  //

  /** 
  if (ReferenceSignal>fData.GetEntries()){
    printf("\n");
    printf("ERROR: ReferenceSignal>nentries %d>%lld \n",ReferenceSignal,fData.GetEntries());
    printf("Reference signal setted to entry with the highest intensity \n");
    printf("\n");  //
  }
  */

  fData.Draw("Entry$",Form("EvtId==%d",ReferenceSignal));
  if (fData.GetSelectedRows()>1)
    cout << " Warning: "<< fData.GetSelectedRows() << " events with the same EvtId " << endl;

  if (fData.GetSelectedRows()==1){ 
    double *yy =  fData.GetV1();
    if (int (yy[0]) > 0) 
      ReferenceSignal = int (yy[0]);
  } else { 
    ReferenceSignal = -100;
  }

  double NeMax = 0.;
  int k_NeMax = 0;
  if (ReferenceSignal<0){
    cout << "Search the signal with the highest Ne ...." << endl;
   for (int k=0; k<fData.GetEntries(); k++) { 
     fData.GetEntry(k);
     if ( SelTriggers(RunId,EvtId,FixVertGain, Ch1,T1, Ch2,T2, Ch3,T3, Ch4,T4, Ne_min, Ne_max, tmin_beam_ped, tmax_beam_ped, tmin_beam_sgn, tmax_beam_sgn, IsDecember2012, SelCh) ){
       GetNe(Ch1, T1, tmin_beam_ped, tmax_beam_ped, tmin_beam_sgn, tmax_beam_sgn, BunchLength, Ne, ped, Ch1_max, Current);
       if (Ne>NeMax){
	 NeMax = Ne;
	 k_NeMax = k;
       } 
     }
   }
   printf("Reference signal:  k=%d  Ne max = %le \n",k_NeMax,NeMax);
   ReferenceSignal = k_NeMax;
  }  

  // time delay defined by alignment
  double deltaT[5] = {0.,0.,0.,0.,0.};

  // tree entry in the list of selected triggers
  int SelReferenceSignal = 0; 

  // save the trace and calculate Ne for the reference signal
  char sdname_alignment[100];
  double NeA = 0.;
  double IA = 0.; // current
  TH1D *h_k_sel_ref = new TH1D("h_k_sel_ref","h_k_sel_ref",1,0.,1.);    // save the index of selected triggers (Entry$ del tree in output)
  TH1D *h_Ne_ref = new TH1D("h_Ne_ref","h_Ne_ref",1,0.,1.);    // save NeA in the output file

  fData.GetEntry(ReferenceSignal);
  GetNe(Ch1, T1, tmin_beam_ped, tmax_beam_ped, tmin_beam_sgn, tmax_beam_sgn, BunchLength, NeA, ped, Ch1_max, IA);
  h_Ne_ref->Fill(0.5,NeA);

  if (Ne_norm>=0. && Ne_norm<10.) Ne_norm = NeA;

  vector<double> TA;
  vector<double> VA;
  vector<double> TbeamA;
  vector<double> VbeamA;
  if (ich_ref==2){
    Put_Pointer_into_Vector(T2,TA);
    Put_Pointer_into_Vector(Ch2,VA);
  }
  if (ich_ref==3){
    Put_Pointer_into_Vector(T3,TA);
    Put_Pointer_into_Vector(Ch3,VA);
  }
  if (ich_ref==4){
    Put_Pointer_into_Vector(T4,TA);
    Put_Pointer_into_Vector(Ch4,VA);
  }

  Put_Pointer_into_Vector(T1,TbeamA);
  VbeamA.clear();
  for (int j=0; j<T1->size(); j++) {
    VbeamA.push_back(Ch1->at(j)+ped);
  }

  printf("\n");
  if (Alignment){ 
    sprintf(sdname_alignment,"Alignment: runid=%d    evtid=%d   Entry$=%d   Ne=%le   Current=%f mA \n",RunId,EvtId,ReferenceSignal,NeA,IA);
    printf("%s\n",sdname_alignment);
    printf("\n");
  }
  int Ncheck = 0;  
  const int max_Ncheck = 100000;
  // define the (bins) range in which we store the traces

  double ttmin = T1->at(0)/ns;
  double tstep = T1->at(1)/ns - T1->at(0)/ns;
  // nota1:  in GetVrms si ha t-deltat > t1  e quindi qui bisogna avere t > t1 + deltat  
  int j1check = ( tmin_plot/ns - ttmin )/tstep; 
  int j2check = ( tmax_plot/ns - ttmin )/tstep; 
  //cout << endl;
  //cout << tmin_plot << " " << tmax_plot << endl;
  //cout << j1check << " " << j2check << endl;
  //cout << endl;
  //int j1check = 4600; 
  //int j2check = 5600; 

  TGraph *gr_vt_check[max_Ncheck];
  TGraph *gr_vt_al_check[max_Ncheck];

  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 



  ////////////////////////////////////////////// frequency cut /////////////////////////////////////////////////////
  MinFreq.clear();  MaxFreq.clear();
  
  if(FreqCutType == 0) {

    // Cut the LINAC Frequency!
    MinFreq.push_back(0.001);  MaxFreq.push_back(3.6);
    MinFreq.push_back(4.9);    MaxFreq.push_back(6.7);
    MinFreq.push_back(7.7);    MaxFreq.push_back(9.5);
    MinFreq.push_back(10.9);   MaxFreq.push_back(12.2);
    MinFreq.push_back(13.8);   MaxFreq.push_back(15.1);
    MinFreq.push_back(16.8);   MaxFreq.push_back(17.5);
    MinFreq.push_back(19.9);   MaxFreq.push_back(19.999);
    
  } else {
    
    
    MinFreq.push_back(0.001);      MaxFreq.push_back(FreqCutMin);
    MinFreq.push_back(FreqCutMax); MaxFreq.push_back(19.999);
    

    /*
    MinFreq.push_back(0.001);  MaxFreq.push_back(3.6);
    MinFreq.push_back(4.9);    MaxFreq.push_back(6.7);
    MinFreq.push_back(7.7);    MaxFreq.push_back(9.5);
    MinFreq.push_back(10.9);   MaxFreq.push_back(12.2);
    MinFreq.push_back(13.8);   MaxFreq.push_back(15.1);
    MinFreq.push_back(16.8);   MaxFreq.push_back(17.5);
    MinFreq.push_back(19.9);   MaxFreq.push_back(19.999);
    */

    /*
    MinFreq.push_back(0.001);  MaxFreq.push_back(4.0);
    MinFreq.push_back(5.0);    MaxFreq.push_back(19.999);
    */
   
    
  }


  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////   
  
  int selfreq = 0;
  double fCalibFactor = GetCalibIntegral(selfreq, MinFreq,  MaxFreq);
  
  selfreq = 1;
  double fCalibFactorCut = GetCalibIntegral(selfreq, MinFreq,  MaxFreq);

  double alpha = 1; //by default
  TF1 *fExpectedMBRVrmsCut = new TF1("fExpectedMBRVrmsCut", fMBRExpectedVRMS, 0.5, 7, 5);
  fExpectedMBRVrmsCut->SetParameter(0, alpha);
  fExpectedMBRVrmsCut->SetParameter(1, Target);
  fExpectedMBRVrmsCut->SetParameter(2, BunchLength);
  fExpectedMBRVrmsCut->SetParameter(3, fCalibFactorCut);
  fExpectedMBRVrmsCut->SetParameter(4, 1);

  TF1 *fExpectedMBRVrms = new TF1("fExpectedMBRVrms", fMBRExpectedVRMS, 0.5, 7, 5);
  fExpectedMBRVrms->SetParameter(0, alpha);
  fExpectedMBRVrms->SetParameter(1, Target);
  fExpectedMBRVrms->SetParameter(2, BunchLength);
  fExpectedMBRVrms->SetParameter(3, fCalibFactor);
  fExpectedMBRVrms->SetParameter(4, 1);


  
  char sdname[100];
  int Naverage_f[Nch][NPeak];
  int Naverage_f_sel[Nch][NPeak];
 
  for (int ich=0; ich<Nch; ich++) { 
    for (int ipeak=0; ipeak<NPeak; ipeak++) { 
      Naverage_f[ich][ipeak] = 0;
      Naverage_f_sel[ich][ipeak] = 0;
     
    }
  }
  // average frequency spectrum
  TH1D *hP_sgn_f[Nch][NPeak], *hP_bkg_f[Nch][NPeak];
  TH1D *hP_sgn_calib_f[Nch][NPeak];
  TH1D *hPdBm_sgn_f[Nch][NPeak], *hPdBm_bkg_f[Nch][NPeak];
  TH1D *hPdBm_sgn_sel_f[Nch][NPeak];

  TH1D *hP_sgn_ref_f[Nch][NPeak], *hP_bkg_ref_f[Nch][NPeak];
  TH1D *hPdBm_sgn_ref_f[Nch][NPeak], *hPdBm_bkg_ref_f[Nch][NPeak];

  int Nbins_sgn_f;
  int Nbins_bkg_f;
  // average time profile
  int Naverage_t = 0;
  TH1D *hBeam_t;
  TH1D *hV_t[Nch];
  TH1D *hV_t_ref[Nch];
  TGraph *hV_t_refCut[Nch];
  TH1D *hP_t[Nch];
  TH1D *hLogPdBm_t[Nch];
  TH1D *hV_vrms_t[Nch];
  TH1D *hP_vrms_t[Nch];
  TH1D *hLogPdBm_vrms_t[Nch];
  double fmin = 0.;
  double fmax = 20.;
  
  /*
  double Power_at_reference[Nch][NPeak];

  for (int ich=0; ich<Nch; ich++) 
    for (int ipeak=0; ipeak<NPeak; ipeak++)
      Power_at_reference[ich][ipeak] = 0;

  */


  //
  // number of events in the data tree
  int nentries = fData.GetEntries();
  //nentries = 1;
  int plot = 1;
  //loop over all events 
  for (int k=0; k<nentries; k++) { 

    // reset tree variables
    InitOutputTree();

    fData.GetEntry(k);
    //printf("\n");

    bool Sel = SelTriggers(RunId,EvtId,FixVertGain, Ch1,T1, Ch2,T2, Ch3,T3, Ch4,T4, Ne_min, Ne_max, tmin_beam_ped, tmax_beam_ped, tmin_beam_sgn, tmax_beam_sgn, IsDecember2012, SelCh);

    if ( Sel ) {

      // Ne bisogna calcolarlo prima di fare l'allineamento in quanto ci serve Ne in AlignSignals
      // comunque Ne non dipende dall'allineamento
      GetNe(Ch1, T1, tmin_beam_ped, tmax_beam_ped, tmin_beam_sgn, tmax_beam_sgn, BunchLength, Ne, ped, Ch1_max, Current);

      // beam changing the time window
      double dt_beam = 20*ns;
      double Current_dummy;
      GetNe(Ch1, T1, tmin_beam_ped, tmax_beam_ped, tmin_beam_sgn-dt_beam, tmax_beam_sgn        , BunchLength, Ne1, ped, Ch1_max, Current_dummy);
      GetNe(Ch1, T1, tmin_beam_ped, tmax_beam_ped, tmin_beam_sgn+dt_beam, tmax_beam_sgn        , BunchLength, Ne2, ped, Ch1_max, Current_dummy);
      GetNe(Ch1, T1, tmin_beam_ped, tmax_beam_ped, tmin_beam_sgn        , tmax_beam_sgn-dt_beam, BunchLength, Ne3, ped, Ch1_max, Current_dummy);
      GetNe(Ch1, T1, tmin_beam_ped, tmax_beam_ped, tmin_beam_sgn        , tmax_beam_sgn+dt_beam, BunchLength, Ne4, ped, Ch1_max, Current_dummy);

      vector<double> *pTA = new vector<double> (TA);
      vector<double> *pVA = new vector<double> (VA);

      vector<double> *pTB = NULL;
      vector<double> *pVB = NULL;

      if (ich_ref==2){
	pTB = T2;
	pVB = Ch2;
      }
      if (ich_ref==3){
	pTB = T3;
	pVB = Ch3;
      }
      if (ich_ref==4){
	pTB = T4;
	pVB = Ch4;
      }

      double dt_ali = 0;
      double chi2 = 0.;
      double NeB = Ne;
      AlignSignals(T_al_1,T_al_2,pVA,pTA,pVB,pTB,NeA,NeB,dt_ali,chi2);
      //dt_ali -= 1.5e-9; //Gonzalo RF 14/10/2013
      for (int ich=1; ich<Nch; ich++) { 
	if (Alignment) {
	  deltaT[ich] = dt_ali;
	} else {
	  deltaT[ich] = 0.;
	}
      }

      // output tree
      if (Alignment) {
	 chi2_al = chi2;
	 dt_al = dt_ali;
      }
      
      if (k%50==0 || k==ReferenceSignal) 
	printf("%4d (%4d/%4d)  runid = %d  evtid = %d   Ne = %4.1le   - All. = %d    dt(align.) = %6.3f ns    chi2(align.) = %5.2f \n",k,Ncheck,nentries,RunId,EvtId,Ne,Alignment,deltaT[1]/ns,chi2_al);    

      if (Ncheck<max_Ncheck){
	gr_vt_check[Ncheck] = new TGraph();
	gr_vt_al_check[Ncheck] = new TGraph();
	for (int j=j1check; j<=j2check; j++) {
	  double x,y;
	  double fac;
	  if (Ne_norm<0.){
	    fac = 1.;
	  } else {
	    fac = Ne_norm/Ne;
	  } 
	  y = pVB->at(j) * fac;
	  x = pTB->at(j);
	  gr_vt_check[Ncheck]->SetPoint(j-j1check,x/ns,y/mV);
	  x = (pTB->at(j) - deltaT[3]);
	  gr_vt_al_check[Ncheck]->SetPoint(j-j1check,x/ns,y/mV);
	}
	if (k==ReferenceSignal) SelReferenceSignal = Ncheck;
	Ncheck++;
      }

      
      // average time profile
      
      Naverage_t++;


      if (Naverage_t==1){
	// if you wish to check the bin time of the oscilloscope
	//stept = T1->at(1)-T1->at(0);
	// to be compared with dt_oscilloscope 
	// time profile histo defined in the same way for all events
	// try to use T4 or T3 or T2 instead of T1
	// the number of bins are the same (see below)
	// however there is a few ps offset (<20 ps) --> see the tree variable TOffset
	int Nbins_t;
	double tmin,tmax;
	Nbins_t = T1->size(); tmin = T1->at(0) - deltaT[1]; tmax = T1->at(Nbins_t-1) - deltaT[1];
	hBeam_t = new TH1D("hBeam_t","hBeam_t",Nbins_t,tmin,tmax); 
      }
	  
      for (unsigned int j=0; j<hBeam_t->GetNbinsX(); j++) {
	// fill graph with beam signal
	// in alternativa il tempo puo' essere definito come
	// x = tmin + j*stept + stept/2.;
	double x,y; 
	x = T1->at(j) - deltaT[1];
	y = Ch1->at(j) + ped;
	// x = tmin + j*stept + stept/2.;
	double fac; 
	if (Ne_norm<0.){
	  fac = 1.;
	} else {
	  fac = Ne_norm/Ne;
	} 
	hBeam_t->Fill(x,y*fac);
      }

      for (int ich=1; ich<Nch; ich++) { 

	if (SelCh[ich]){

	  vector<double> *Time = NULL; 
	  vector<double> *Ch = NULL;

	  if (ich==2){ 
	    Time = T2; 
	    Ch = Ch2; 
	  }
	  if (ich==3){ 
	    Time = T3; 
	    Ch = Ch3; 
	  }
	  if (ich==4){ 
	    Time = T4; 
	    Ch = Ch4; 
	  }

	  // trace with the signal time selection
	  if (Naverage_t==1){
	    gr_vt[ich] = new TGraph();
	    gr_vt_sel0[ich] = new TGraph();
	    gr_vt_sel1[ich] = new TGraph();
	    gr_vt_sel2[ich] = new TGraph();
	  }
	  if (k==ReferenceSignal) {
	    int jj = 0;
	    int jj0 = 0;
	    int jj1 = 0;
	    int jj2 = 0;
	    for (int j=j1check; j<=j2check; j++) {
	      double x,y;
	      x = Time->at(j);
	      double fac; 
	      if (Ne_norm<0.){
		fac = 1.;
	      } else {
		fac = Ne_norm/Ne;
	      } 
	      y = Ch->at(j)*fac;
	      gr_vt[ich]->SetPoint(jj,x/ns,y/mV); 
	      jj++;
	      if ( x>tmin_sgn[ich] && x<tmax_sgn[ich] ){
		gr_vt_sel0[ich]->SetPoint(jj0,x/ns,y/mV); 
		jj0++;
	      }
	      if ( x>t1min_sgn[ich] && x<t1max_sgn[ich] ){
		gr_vt_sel1[ich]->SetPoint(jj1,x/ns,y/mV); 
		jj1++;
	      }
	      if ( x>t2min_sgn[ich] && x<t2max_sgn[ich] ){
		gr_vt_sel2[ich]->SetPoint(jj2,x/ns,y/mV); 
		jj2++;	      }
	    }
	  }

	  if (Naverage_t==1){
	    char dummy[150];
	    int Nbins_t;
	    double tmin,tmax;
	    Nbins_t = Time->size(); tmin = (Time->at(0) - deltaT[ich])/ns; tmax = (Time->at(Nbins_t-1) - deltaT[ich])/ns;
	    strcpy (sdname,"hV_sgn"); sprintf(dummy,"%d",ich); strcat (sdname,dummy); strcat (sdname,"_t"); 
	    hV_t[ich]= new TH1D(sdname,sdname,Nbins_t,tmin,tmax);

	    strcpy (sdname,"hV_ref_sgn"); sprintf(dummy,"%d",ich); strcat (sdname,dummy); strcat (sdname,"_t"); 
	    hV_t_ref[ich]= new TH1D(sdname,sdname,Nbins_t,tmin,tmax);
	    
	    //strcpy (sdname,"hV_refCut_sgn"); sprintf(dummy,"%d",ich); strcat (sdname,dummy); strcat (sdname,"_t"); 
	    hV_t_refCut[ich]= new TGraph();

	    strcpy (sdname,"hP_sgn"); sprintf(dummy,"%d",ich); strcat (sdname,dummy); strcat (sdname,"_t"); 
	    hP_t[ich]= new TH1D(sdname,sdname,Nbins_t,tmin,tmax);
	    strcpy (sdname,"hLogPdBm_sgn"); sprintf(dummy,"%d",ich); strcat (sdname,dummy); strcat (sdname,"_t"); 
	    hLogPdBm_t[ich]= new TH1D(sdname,sdname,Nbins_t,tmin,tmax);

	    strcpy (sdname,"hV_sgn"); sprintf(dummy,"%d",ich); strcat (sdname,dummy); strcat (sdname,"_vrms_t"); 
	    hV_vrms_t[ich]= new TH1D(sdname,sdname,Nbins_t,tmin,tmax);
	    strcpy (sdname,"hP_sgn"); sprintf(dummy,"%d",ich); strcat (sdname,dummy); strcat (sdname,"_vrms_t"); 
	    hP_vrms_t[ich]= new TH1D(sdname,sdname,Nbins_t,tmin,tmax);
	    strcpy (sdname,"hLogPdBm_sgn"); sprintf(dummy,"%d",ich); strcat (sdname,dummy); strcat (sdname,"_vrms_t"); 
	    hLogPdBm_vrms_t[ich]= new TH1D(sdname,sdname,Nbins_t,tmin,tmax);
	  }
	
	 
	  for (unsigned int j=0; j<hV_t[ich]->GetNbinsX(); j++) {
	    double fac; 
	    if (Ne_norm<0.){
	      fac = 1.;
	    } else {
	      fac = Ne_norm/Ne;
	    } 
	    double x,y;
	    x = (Time->at(j) - deltaT[ich])/ns;
	    y = Ch->at(j);
	    hV_t[ich]->Fill(x,y/mV*fac);
	    if (k==ReferenceSignal)  
	      hV_t_ref[ich]->Fill(x,y/mV*fac);
	    
	    y = Ch->at(j) * Ch->at(j) / R;
	    hP_t[ich]->Fill(x,y*fac*fac);

	    double dt = 1.; 
	    double tmin = x - dt/2.;
	    double tmax = x + dt/2.;
	    y = GetVrms_bin(Ch, Time, deltaT[ich], tmin*ns, tmax*ns);
	    hV_vrms_t[ich]->Fill(x,y/mV*fac);
	    hP_vrms_t[ich]->Fill(x,y*y/R*fac*fac);

	  }

	}

      }




      for (int ich=1; ich<Nch; ich++) { 

	if (SelCh[ich]){

	  vector<double> *Time = NULL; 
	  vector<double> *Ch = NULL;

	  if (ich==2){ 
	    Time = T2; 
	    Ch = Ch2; 
	  }
	  if (ich==3){ 
	    Time = T3; 
	    Ch = Ch3; 
	  }
	  if (ich==4){ 
	    Time = T4; 
	    Ch = Ch4; 
	  }

	  for (int ipeak=0; ipeak<NPeak; ipeak++) { 

	    if (SelPeak[ipeak]){

	      // time window definition
	      double tmin,tmax;
	      if (ipeak==0){
		tmin = tmin_sgn[ich];
		tmax = tmax_sgn[ich];
	      }
	      if (ipeak==1){
		tmin = t1min_sgn[ich];
		tmax = t1max_sgn[ich];
	      }
	      if (ipeak==2){
		tmin = t2min_sgn[ich];
		tmax = t2max_sgn[ich];
	      }

	      //printf("%d: peak=%d  t [%le %le]  [%le,%le] \n",Naverage_f[ich][ipeak],ipeak,tmin,tmax,Ti1,Tf1, t1min_sgn[ich], t1max_sgn[5]);

	      Fsgn.clear(); Psgn.clear(); VCutsgn.clear(); TCutsgn.clear(); PCutsgn.clear();
	      Fbkg.clear(); Pbkg.clear(); VCutbkg.clear(); TCutbkg.clear(); PCutbkg.clear();
	      TransformFFT(Ch, Time, deltaT[ich], tmin    , tmax    , MinFreq, MaxFreq, Fsgn, Psgn, VCutsgn, TCutsgn, PCutsgn);
	      TransformFFT(Ch, Time, deltaT[ich], tmin_bkg, tmax_bkg, MinFreq, MaxFreq, Fbkg, Pbkg, VCutbkg, TCutbkg, PCutbkg);
	     
	      // output tree
	      

// simplify the code Valerio 13/04/2014

	      double vrms_sgn = GetVrms(Ch, Time, deltaT[ich], tmin,     tmax);
	      double vrms_bkg = GetVrms(Ch, Time, deltaT[ich], tmin_bkg, tmax_bkg     );
	      double vrms_sgn_fsel1 =  GetVrms1(VCutsgn, TCutsgn, deltaT[ich], tmin,     tmax);
	      double vrms_bkg_fsel1 =  GetVrms1(VCutbkg, TCutbkg, deltaT[ich], tmin_bkg, tmax_bkg     );
	      fExpectedMBRVrms->SetParameter(0,alpha=1);
	      fExpectedMBRVrmsCut->SetParameter(0,alpha=1);
	      double vrms_MBR1 = fExpectedMBRVrms->Eval(Ne/1e9);
	      double vrms_MBR1_fsel1 = fExpectedMBRVrmsCut->Eval(Ne/1e9);
	      fExpectedMBRVrmsCut->SetParameter(0,alpha=2);
	      fExpectedMBRVrms->SetParameter(0,alpha=2);
	      double vrms_MBR2 = fExpectedMBRVrms->Eval(Ne/1e9);
	      double vrms_MBR2_fsel2 = fExpectedMBRVrmsCut->Eval(Ne/1e9);

	      if (ich==1){
		vrms_sgn_ch1[ipeak]       = vrms_sgn;
		vrms_bkg_ch1[ipeak]       = vrms_bkg;
		vrms_sgn_ch1_fsel1[ipeak] = vrms_sgn_fsel1;
		vrms_bkg_ch1_fsel1[ipeak] = vrms_bkg_fsel1;
		vrms_dedx_exp1_sgn_ch1[ipeak]       = vrms_MBR1;    
		vrms_dedx_exp1_sgn_ch1_fsel1[ipeak] = vrms_MBR1_fsel1; 
		vrms_dedx_exp2_sgn_ch1[ipeak]       = vrms_MBR2;
		vrms_dedx_exp2_sgn_ch1_fsel1[ipeak] = vrms_MBR2_fsel2;
	      } 
	      if (ich==2){
		vrms_sgn_ch2[ipeak]       = vrms_sgn;
		vrms_bkg_ch2[ipeak]       = vrms_bkg;
		vrms_sgn_ch2_fsel1[ipeak] = vrms_sgn_fsel1;
		vrms_bkg_ch2_fsel1[ipeak] = vrms_bkg_fsel1;
		vrms_dedx_exp1_sgn_ch2[ipeak]       = vrms_MBR1;    
		vrms_dedx_exp1_sgn_ch2_fsel1[ipeak] = vrms_MBR1_fsel1; 
		vrms_dedx_exp2_sgn_ch2[ipeak]       = vrms_MBR2;
		vrms_dedx_exp2_sgn_ch2_fsel1[ipeak] = vrms_MBR2_fsel2;
	      } 
	      if (ich==3){
		vrms_sgn_ch3[ipeak]       = vrms_sgn;
		vrms_bkg_ch3[ipeak]       = vrms_bkg;
		vrms_sgn_ch3_fsel1[ipeak] = vrms_sgn_fsel1;
		vrms_bkg_ch3_fsel1[ipeak] = vrms_bkg_fsel1;
		vrms_dedx_exp1_sgn_ch3[ipeak]       = vrms_MBR1;    
		vrms_dedx_exp1_sgn_ch3_fsel1[ipeak] = vrms_MBR1_fsel1; 
		vrms_dedx_exp2_sgn_ch3[ipeak]       = vrms_MBR2;
		vrms_dedx_exp2_sgn_ch3_fsel1[ipeak] = vrms_MBR2_fsel2;
	      } 
	      if (ich==4){
		vrms_sgn_ch4[ipeak]       = vrms_sgn;
		vrms_bkg_ch4[ipeak]       = vrms_bkg;
		vrms_sgn_ch4_fsel1[ipeak] = vrms_sgn_fsel1;
		vrms_bkg_ch4_fsel1[ipeak] = vrms_bkg_fsel1;
		//Power_at_reference[ich][ipeak] = ( pow(vrms_sgn_fsel1,2) - pow(vrms_bkg_fsel1,2) ) / 50;
		vrms_dedx_exp1_sgn_ch4[ipeak]       = vrms_MBR1;    
		vrms_dedx_exp1_sgn_ch4_fsel1[ipeak] = vrms_MBR1_fsel1; 
		vrms_dedx_exp2_sgn_ch4[ipeak]       = vrms_MBR2;
		vrms_dedx_exp2_sgn_ch4_fsel1[ipeak] = vrms_MBR2_fsel2;
	      } 

	      // che cosa e'?
	      // to save the Aeff ... plots ?
	      if (plot == 1) {
		ExpectedVrms_MBRdEdX(2.5e9, 1, 3, Target, BunchLength, 0, Time, deltaT[ich],  tmin,  tmax,  MinFreq, MaxFreq);
		plot = 100;
	      }

	     	      
	      if (k==ReferenceSignal) {
		double fac; 
		if (Ne_norm<0.){
		  fac = 1.;
		} else {
		  fac = Ne_norm/Ne;
		} 
		for (unsigned int j=0; j<TCutsgn.size(); ++j) {
		  double x = (TCutsgn.at(j) - deltaT[ich])/ns;
		  double y = VCutsgn.at(j);
		  hV_t_refCut[ich]->SetPoint(j,x,y/mV*fac);
		}
	      }
	      
	    // Average frequency spectrum

	      Naverage_f[ich][ipeak]++;

	      if (Naverage_f[ich][ipeak]==1){

		Nbins_sgn_f = Fsgn.size();
		Nbins_bkg_f = Fbkg.size();
		//printf("Ch%d  sgn %d (%le %le)   bkg %d (%le %le)\n",ich,Nbins_sgn_f[ich],tmin_sgn[ich],tmax_sgn[ich],Nbins_bkg_f,tmin_bkg,tmax_bkg);
		char dummy[150];

		strcpy (sdname,"hP_sgn"); sprintf(dummy,"%d_ref_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_f"); 
		hP_sgn_ref_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_sgn_f,fmin,fmax);
		strcpy (sdname,"hP_bkg"); sprintf(dummy,"%d_ref_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_f");
		hP_bkg_ref_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_bkg_f,fmin,fmax);

		strcpy (sdname,"hPdBm_sgn"); sprintf(dummy,"%d_ref_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_f");
		hPdBm_sgn_ref_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_sgn_f,fmin,fmax);
		strcpy (sdname,"hPdBm_bkg"); sprintf(dummy,"%d_ref_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_f");
		hPdBm_bkg_ref_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_bkg_f,fmin,fmax);

		strcpy (sdname,"hP_sgn"); sprintf(dummy,"%d_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_f"); 
		hP_sgn_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_sgn_f,fmin,fmax);
		strcpy (sdname,"hP_bkg"); sprintf(dummy,"%d_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_f");
		hP_bkg_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_bkg_f,fmin,fmax);
		
		strcpy (sdname,"hPdBm_sgn"); sprintf(dummy,"%d_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_f");
		hPdBm_sgn_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_sgn_f,fmin,fmax);
		strcpy (sdname,"hPdBm_bkg"); sprintf(dummy,"%d_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_f");
		hPdBm_bkg_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_bkg_f,fmin,fmax);
		
		strcpy (sdname,"hP_sgn"); sprintf(dummy,"%d_calib_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_f"); 
		hP_sgn_calib_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_sgn_f,fmin,fmax);
		strcpy (sdname,"hPdBm_sgn"); sprintf(dummy,"%d_peak%d",ich,ipeak); strcat (sdname,dummy); strcat (sdname,"_sel_f");
		hPdBm_sgn_sel_f[ich][ipeak] = new TH1D(sdname,sdname,Nbins_sgn_f,fmin,fmax);

	      }

	      Nbins_sgn_f = Fsgn.size();
	      Nbins_bkg_f = Fbkg.size();

	      if ( Nbins_sgn_f != hP_sgn_f[ich][ipeak]->GetNbinsX()|| 
		   Nbins_bkg_f != hP_bkg_f[ich][ipeak]->GetNbinsX()         ){
		printf("WARNING Ch%d sgn %ld %d %d\n",ich,Psgn.size(),Nbins_sgn_f,hP_sgn_f[ich][ipeak]->GetNbinsX());
		printf("WARNING Ch%d bkg %ld %d %d\n",ich,Pbkg.size(),Nbins_bkg_f,hP_bkg_f[ich][ipeak]->GetNbinsX());
		
	      } else {
	  
		Naverage_f_sel[ich][ipeak]++;


		// per forzare la normalizzazione corretta
		// le due opzioni sono simili (al %) ma non esattamente uguali
		// Psgn.at(j)/R * / hP_sgn_f[ich][ipeak]->GetNbinsX();
		// Psgn.at(j) * fac_sgn_true
		double ytot;
		
		ytot = 0;
		for (int j=0; j<Nbins_sgn_f; j++) 
		  ytot = ytot + Psgn.at(j);
		double fac_sgn_true = pow(vrms_sgn,2)/R/ytot;
		
		ytot = 0;
		for (int j=0; j<Nbins_bkg_f; j++) 
		  ytot = ytot + Pbkg.at(j);
		double fac_bkg_true = pow(vrms_bkg,2)/R/ytot;



 		for (int j=0; j<Nbins_sgn_f; j++) {
		  double fac; 
		  if (Ne_norm<0.){
		    fac = 1.;
		  } else {
		    fac = Ne_norm/Ne;
		    // dovrebbe essere 
		    // in realta' la potenza andrebbe riscalata nel modo seguente:  P(Ne') = (P(Ne)-P_bkg) * (Ne'/Ne)^2 + P_bkg
		    fac = pow(Ne_norm/Ne,2);

		  } 

		  //
		  fac = 1;
		  //

		  double step = (fmax-fmin)/( double(Nbins_sgn_f) );
		  double nu = fmin + j*step + step/2.;
		  double x,y; 
		  x = nu;
		  //y = Psgn.at(j)/R * fac;
		  y = Psgn.at(j)/R * fac / hP_sgn_f[ich][ipeak]->GetNbinsX();
		  //y = Psgn.at(j) * fac_sgn_true;
		  hP_sgn_f[ich][ipeak]->Fill(x,y);
		  if (k==ReferenceSignal) 
		    hP_sgn_ref_f[ich][ipeak]->Fill(x,y);
		  
		}
	    
		for (int j=0; j<Nbins_bkg_f; j++) {
		  double fac; 
		  if (Ne_norm<0.){
		    fac = 1.;
		  } else {
		    fac = Ne_norm/Ne;
		  } 
		  // il fondo non puo' dipendere da Ne
		  fac = 1;
		  //
		  double step = (fmax-fmin)/( double(Nbins_bkg_f) );
		  double nu = fmin + j*step + step/2.;
		  double x,y; 
		  x = nu;
		  //y = Pbkg.at(j)/R * fac;
		  y = Pbkg.at(j)/R * fac / hP_sgn_f[ich][ipeak]->GetNbinsX();
		  hP_bkg_f[ich][ipeak]->Fill(x,y);
		  if (k==ReferenceSignal) 
		    hP_bkg_ref_f[ich][ipeak]->Fill(x,y);
		}
	      }


	    } // peak selection

	  } // peak loop

	} // channel selection

      } // channel loop
     
      // output tree
      runid = RunId;
      evtid = EvtId;
      bunch_length = BunchLength;
      target = Target;
      //Ne = Ne;      
      //Ne1 = Ne1;      
      //Ne2 = Ne2;      
      //Ne3 = Ne3;      
      //Ne4 = Ne4;      
      //Current = Current
      fTreeOut->Fill();

     } // trigger selection

  } // loop over triggers

  // Normalise the average time profile
  hBeam_t->Scale(1./Naverage_t);
  for (int ich=1; ich<Nch; ich++) { 
    if (SelCh[ich]){
      hV_t[ich]->Scale(1./Naverage_t);
      hP_t[ich]->Scale(1./Naverage_t);
      for (int l=1; l<=hP_t[ich]->GetNbinsX(); l++) { 
	double y = hP_t[ich]->GetBinContent(l);
	//if (y>0) y = log10( y);
	if (y>0) y = 10.*log10( y*1e3 );
	hLogPdBm_t[ich]->SetBinContent(l,y);
      }
      hV_vrms_t[ich]->Scale(1./Naverage_t);
      hP_vrms_t[ich]->Scale(1./Naverage_t);
      for (int l=1; l<=hP_vrms_t[ich]->GetNbinsX(); l++) { 
	double y = hP_vrms_t[ich]->GetBinContent(l);
	//if (y>0) y = log10( y);
	if (y>0) y = 10.*log10( y*1e3 );
	hLogPdBm_vrms_t[ich]->SetBinContent(l,y);
      }
    }
  }
  

  // average signal after a frequency cut

  TGraph *grV_sgn_Cut_t[Nch], *grPdBm_sgn_Cut_f[Nch], *grPdBm_bkg_Cut_f[Nch];

  //MinFreq.clear();  MaxFreq.clear();
  //MinFreq.push_back(Fmin_av);  MaxFreq.push_back(Fmax_av);  


  for (int ich=1; ich<Nch; ich++) 
    if (SelCh[ich]) {
      
      vector<double>  fTimeAvg;
      vector<double>  fVoltAvg;
      fTimeAvg.clear(), fVoltAvg.clear();
      for (int l=1; l<=hV_t[ich]->GetNbinsX(); l++) { 
	double tmp_time = hV_t[ich]->GetBinCenter(l) * ns; 
	double tmp_volt = hV_t[ich]->GetBinContent(l) * mV;
	fTimeAvg.push_back(tmp_time);
	fVoltAvg.push_back(tmp_volt);
      }

      Fsgn.clear(); Psgn.clear(); VCutsgn.clear(); TCutsgn.clear(); PCutsgn.clear();
      Fbkg.clear(); Pbkg.clear(); VCutbkg.clear(); TCutbkg.clear(); PCutbkg.clear();
      TransformFFT(&fVoltAvg, &fTimeAvg, 0, Ti_all, Tf_all,     MinFreq, MaxFreq, Fsgn, Psgn, VCutsgn, TCutsgn, PCutsgn);
      TransformFFT(&fVoltAvg, &fTimeAvg, 0, tmin_bkg, tmax_bkg, MinFreq, MaxFreq, Fbkg, Pbkg, VCutbkg, TCutbkg, PCutbkg);
      
      grV_sgn_Cut_t[ich] = new TGraph(); 
      grPdBm_sgn_Cut_f[ich] = new TGraph(); 
      grPdBm_bkg_Cut_f[ich] = new TGraph(); 
      
      for (int t=0;t<TCutsgn.size();++t)
	grV_sgn_Cut_t[ich]->SetPoint(t,TCutsgn[t]/ns,VCutsgn[t]/mV);
      
      for (int f=0;f<Fsgn.size();++f) {
	double y = Psgn[f]/R;
	if (y>0) y = 10 * log10(y*1e3);
	grPdBm_sgn_Cut_f[ich]->SetPoint(f,Fsgn[f]/GHz,y);
      }
      for (int f=0;f<Fbkg.size();++f) {
	double y = Pbkg[f]/R;
	if (y>0) y = 10 * log10(y*1e3);
	grPdBm_bkg_Cut_f[ich]->SetPoint(f,Fbkg[f]/GHz,y);
      }      
      
    }
  

  // Normalise the average frequency spectrum
  for (int ich=1; ich<Nch; ich++) { 
    if (SelCh[ich]){
      for (int ipeak=0; ipeak<NPeak; ipeak++) { 
	if (SelPeak[ipeak]){
	  hP_sgn_f[ich][ipeak]->Scale(1./Naverage_f_sel[ich][ipeak]);
	 
	  for (int l=1; l<=hP_sgn_f[ich][ipeak]->GetNbinsX(); l++) { 
	    double y = hP_sgn_f[ich][ipeak]->GetBinContent(l);
	    if (y>0) y = 10.*log10( y*1e3 );
	    hPdBm_sgn_f[ich][ipeak]->SetBinContent(l,y);
	  }

	  hP_bkg_f[ich][ipeak]->Scale(1./Naverage_f_sel[ich][ipeak]);
	  for (int l=1; l<=hP_bkg_f[ich][ipeak]->GetNbinsX(); l++) { 
	    double y = hP_bkg_f[ich][ipeak]->GetBinContent(l);
	    if (y>0) y = 10.*log10( y*1e3 );
	    hPdBm_bkg_f[ich][ipeak]->SetBinContent(l,y);
	  }

	  for (int l=1; l<=hP_sgn_ref_f[ich][ipeak]->GetNbinsX(); l++) { 
	    double y = hP_sgn_ref_f[ich][ipeak]->GetBinContent(l);
	    if (y>0) y = 10.*log10( y*1e3 );
	    hPdBm_sgn_ref_f[ich][ipeak]->SetBinContent(l,y);
	  }

	  for (int l=1; l<=hP_bkg_ref_f[ich][ipeak]->GetNbinsX(); l++) { 
	    double y = hP_bkg_ref_f[ich][ipeak]->GetBinContent(l);
	    if (y>0) y = 10.*log10( y*1e3 );
	    hPdBm_bkg_ref_f[ich][ipeak]->SetBinContent(l,y);
	  }


	  // spectrum in the selected frequency range
	  for (int j=0; j<hPdBm_sgn_f[ich][ipeak]->GetNbinsX(); j++) {
	    double x = hPdBm_sgn_f[ich][ipeak]->GetBinCenter(j+1);
	    double y = hPdBm_sgn_f[ich][ipeak]->GetBinContent(j+1);
	    bool Sel_f = true;
	    for (int jf=0; jf<MinFreq.size(); jf++) {
	      if (x>MinFreq.at(jf) && x<=MaxFreq.at(jf))
		Sel_f = false;
	    }
	    if (Sel_f) {
	      hPdBm_sgn_sel_f[ich][ipeak]->SetBinContent(j+1,y);
	    } else { 
	      hPdBm_sgn_sel_f[ich][ipeak]->SetBinContent(j+1,-200.);
	    }
	  }
	}
      }
    }
  }



  // write the output

  cout << endl;

  // Alignment 
  TGraph *gr_vt_ref, *gr_vt_ref_dt_al;
  TCanvas *c_vt, *c_vt_al;
  /*
  for (int l=0; l<Ncheck; l++) {
    sprintf(sdname,"gr_vt_check_%d",l); 
    gr_vt_check[l]->Write(sdname);
    sprintf(sdname,"gr_vt_al_check_%d",l); 
    gr_vt_al_check[l]->Write(sdname);
  }
  */

  // store separately the reference signal 
  h_Ne_ref->Write();
  h_k_sel_ref->Fill(0.5,SelReferenceSignal);
  h_k_sel_ref->Write();
  gr_vt_ref = new TGraph();
  gr_vt_ref = gr_vt_check[SelReferenceSignal];
  // SelReferenceSignal salvato in h_k_sel_ref
  //sprintf(sdname,"gr_vt_ref_%d",SelReferenceSignal); 
  sprintf(sdname,"gr_vt_ref");
  gr_vt_ref->Write(sdname);
  // store the reference signal in the time range selected to make the alignment
  gr_vt_ref_dt_al = new TGraph();
  int N_ref_dt_al = 0; 
  for (int j=0; j<gr_vt_ref->GetN(); j++) {
    double x,y;
    gr_vt_ref->GetPoint(j,x,y);
    if (x>(T_al_1/ns) && x<(T_al_2/ns)){
      gr_vt_ref_dt_al->SetPoint(N_ref_dt_al,x,y);
      N_ref_dt_al++;
    }
  }
  //sprintf(sdname,"gr_vt_ref_dt_al_%d",SelReferenceSignal); 
  sprintf(sdname,"gr_vt_ref_dt_al");
  gr_vt_ref_dt_al->SetLineColor(4);
  gr_vt_ref_dt_al->Write(sdname);
  // overimpose all signals
  // not aligned
  c_vt = new TCanvas("c_vt","c_vt",950,850);
  for (int l=0; l<Ncheck; l++) {
    if (l==0) gr_vt_check[l]->Draw("AL");
    if (l>0)  gr_vt_check[l]->Draw("LSAME");
  }
  c_vt->Write();
  // aligned
  c_vt_al = new TCanvas("c_vt_al","c_vt_al",950,850);
  for (int l=0; l<Ncheck; l++) {
    if (l==0) gr_vt_al_check[l]->Draw("AL");
    if (l>0) gr_vt_al_check[l]->Draw("LSAME");
  }
  c_vt_al->Write();


  // trace with the signal in the selected time window
  for (int ich=1; ich<Nch; ich++) { 
    if (SelCh[ich]){
      sprintf(sdname,"gr_vt_ch%d",ich); 
      gr_vt[ich]->Write(sdname);
      sprintf(sdname,"gr_vt_sel0_ch%d",ich); 
      gr_vt_sel0[ich]->Write(sdname);
      sprintf(sdname,"gr_vt_sel1_ch%d",ich); 
      gr_vt_sel1[ich]->Write(sdname);
      sprintf(sdname,"gr_vt_sel2_ch%d",ich); 
      gr_vt_sel2[ich]->Write(sdname);
    }
  }

  // beam plots
  hBeam_t->Write();
  TGraph *gr_AvBeam = new TGraph();
  int N_gr_AvBeam = 0;
  TGraph *gr_AvBeam_sgn = new TGraph();
  int N_gr_AvBeam_sgn = 0;
  for (int j=1; j<hBeam_t->GetNbinsX(); j++) { 
    double x,y;
    x = hBeam_t->GetBinCenter(j);
    y = hBeam_t->GetBinContent(j);
    gr_AvBeam->SetPoint(N_gr_AvBeam,x/ns,y);
    N_gr_AvBeam++;
    if (x>tmin_beam_sgn && x<tmax_beam_sgn) {
      gr_AvBeam_sgn->SetPoint(N_gr_AvBeam_sgn,x/ns,y);
      N_gr_AvBeam_sgn++;
    }
  }
  gr_AvBeam->Write("gr_AvBeam");
  gr_AvBeam_sgn->Write("gr_AvBeam_sgn");
  TGraph *gr_Beam_ref = new TGraph();
  for (int j=0; j<TbeamA.size(); j++) {
    gr_Beam_ref->SetPoint(j,TbeamA.at(j)/ns,VbeamA.at(j));
  }
  gr_Beam_ref->Write("gr_Beam_ref");


  for (int ich=1; ich<Nch; ich++) { 
    if (SelCh[ich]){
      hV_t[ich]->Write();
      hV_t_ref[ich]->Write();
      char dummy[150];
      strcpy (sdname,"hV_refCut_sgn"); sprintf(dummy,"%d",ich); strcat (sdname,dummy); strcat (sdname,"_t"); 
      hV_t_refCut[ich]->Write(sdname);
      hP_t[ich]->Write();
      hLogPdBm_t[ich]->Write();
      hV_vrms_t[ich]->Write();
      hP_vrms_t[ich]->Write();
      hLogPdBm_vrms_t[ich]->Write();

      grV_sgn_Cut_t[ich]->Write(Form("grV_sgn%d_Cut_t",ich));
      grPdBm_sgn_Cut_f[ich]->Write(Form("grPdBm_sgn%d_Cut_f",ich));
      grPdBm_bkg_Cut_f[ich]->Write(Form("grPdBm_bkg%d_Cut_f",ich));

    }
  }

  for (int ich=1; ich<Nch; ich++) {
    if (SelCh[ich]){
      for (int ipeak=0; ipeak<NPeak; ipeak++) { 
	if (SelPeak[ipeak]){
	  hP_sgn_ref_f[ich][ipeak]->Write();
	  hP_bkg_ref_f[ich][ipeak]->Write();
	  hPdBm_sgn_ref_f[ich][ipeak]->Write();
	  hPdBm_bkg_ref_f[ich][ipeak]->Write();
	}
      }
    }
  }



  for (int ich=1; ich<Nch; ich++) {
    if (SelCh[ich]){
      for (int ipeak=0; ipeak<NPeak; ipeak++) { 
	if (SelPeak[ipeak]){
	  hP_sgn_f[ich][ipeak]->Write();
	  hP_bkg_f[ich][ipeak]->Write();

	  
	  for (int j=0; j<hP_sgn_f[ich][ipeak]->GetNbinsX(); j++) {
	      
	      double nu    = hP_sgn_f[ich][ipeak]->GetBinCenter(j+1);
	      double y_sgn = hP_sgn_f[ich][ipeak]->GetBinContent(j+1);
	      
	      double y_bkg = hP_bkg_f[ich][ipeak]->Interpolate(nu);
	      
	      //if(ich==3)
	      //y_bkg = hP_bkg_f[ich][ipeak]->GetBinContent(j+1);
	      
	      
	      double Gain = tgAntenna->Eval(nu*GHz);
	      double y_calib = (y_sgn - y_bkg)/Calib_Dec2012(nu);
	      
	      hP_sgn_calib_f[ich][ipeak]->SetBinContent(j+1,y_calib);
	    }

	  hP_sgn_calib_f[ich][ipeak]->Write();

	  hPdBm_sgn_f[ich][ipeak]->Write();
	  hPdBm_bkg_f[ich][ipeak]->Write();
	  hPdBm_sgn_sel_f[ich][ipeak]->Write();
	}
      }
    }
  }


  
  TCanvas *cPower = new TCanvas("cPower","",640,480);

  double bound_ne = 0;
  for(int b=1;b<1000;++b) {
    
    bound_ne = double(b)/1000.0;
    fTreeOut->Draw("Ne/1e9", Form("Ne/1e9 > %f && Ne/1e9 < %f",(1-bound_ne) * NeA/1e9, (1+bound_ne) * NeA/1e9),"");
    
    TH1F * thNe =  (TH1F*)cPower->FindObject("htemp");
    if(thNe->GetEntries() > 100)
      break;
    

  }
  

  TGraph *tgMBRFluxLimitFull[NPeak];
  TH1F * thPowerFull[NPeak];
  for(int ich=1;ich<Nch;++ich)
    if (SelCh[ich])
      for(int ipeaks=0;ipeaks<NPeak;++ipeaks) {
	
	cout << "Full bandwidth ..." << endl;

	fTreeOut->Draw(Form("(pow(vrms_sgn_ch3[%d],2) - pow(vrms_bkg_ch3[%d],2))/50 * 1e6",ipeaks,ipeaks),
		       Form("Ne/1e9 > %f && Ne/1e9 < %f",(1-bound_ne) * NeA/1e9, (1+bound_ne) * NeA/1e9),"");

	thPowerFull[ipeaks] = (TH1F*)cPower->FindObject("htemp");	

	cout << "channel: " << ich << "\t window: " << ipeaks << "\t Ne: " << NeA/1e9 << "\t P_full " << thPowerFull[ipeaks]->GetMean() << endl;

	tgMBRFluxLimitFull[ipeaks] = new TGraph();

	//alpha value from [1,2]
	int NN = 20;
	double alpha_min = 1, alpha_max=2;
	for(int i=0;i<NN;++i)
	  {

	    double slope = alpha_min + (alpha_max - alpha_min) * i /(NN-1);
	    fExpectedMBRVrms->FixParameter(0, slope);
	    double exp_Vrms = fExpectedMBRVrms->Eval(NeA/1e9);
	    double p_exp = pow(exp_Vrms,2)/50. * 1e6;
	    double I_limit = thPowerFull[ipeaks]->GetMean()/p_exp;
	    tgMBRFluxLimitFull[ipeaks]->SetPoint(i,slope,I_limit);



	  }

	if(ipeaks == 0) {
		
	  tgMBRFluxLimitFull[ipeaks]->SetLineColor(1);
	  tgMBRFluxLimitFull[ipeaks]->SetFillColor(1);

	} else {

	  tgMBRFluxLimitFull[ipeaks]->SetLineColor(ipeaks * 2);
	  tgMBRFluxLimitFull[ipeaks]->SetFillColor(ipeaks * 2);

	}


	tgMBRFluxLimitFull[ipeaks]->SetLineWidth(2);
	tgMBRFluxLimitFull[ipeaks]->SetLineStyle(2);
	tgMBRFluxLimitFull[ipeaks]->SetTitle(Form("MBR Intensity at Ne/1e9 = %f",NeA/1e9));
	tgMBRFluxLimitFull[ipeaks]->GetYaxis()->SetTitle("I/I_{0} [I_{0}=4#times 10^{-16} W/m^{2}/Hz]");
	tgMBRFluxLimitFull[ipeaks]->GetXaxis()->SetTitle("alpha");

	tgMBRFluxLimitFull[ipeaks]->Write(Form("tgMBRFluxLimitFull_sgn%d_peak%d",ich,ipeaks));
	thPowerFull[ipeaks]->Write(Form("thPowerFull_sgn%d_peak%d",ich,ipeaks));


      }


  TGraph *tgMBRFluxLimitCut[NPeak];
  TH1F * thPowerCut[NPeak];
  for(int ich=1;ich<Nch;++ich)
    if (SelCh[ich])
      for(int ipeaks=0;ipeaks<NPeak;++ipeaks) {
	
	cout << "Cut bandwidth ..." << endl;
	
	fTreeOut->Draw(Form("(pow(vrms_sgn_ch3_fsel1[%d],2) - pow(vrms_bkg_ch3_fsel1[%d],2))/50 * 1e6",ipeaks,ipeaks),
		       Form("Ne/1e9 > %f && Ne/1e9 < %f",(1-bound_ne) * NeA/1e9, (1+bound_ne) * NeA/1e9),"");
	
	thPowerCut[ipeaks] = (TH1F*)cPower->FindObject("htemp");

	cout << "channel: " << ich << "\t window: " << ipeaks << "\t Ne: " << NeA/1e9 << "\t P_cut: " << thPowerCut[ipeaks]->GetMean() << endl;

	tgMBRFluxLimitCut[ipeaks] = new TGraph();
	
	//alpha value from [1,2]
	int NN = 20;
	double alpha_min = 1, alpha_max=2;
	for(int i=0;i<NN;++i)
	  {

	    
	    double slope = alpha_min + (alpha_max - alpha_min) * i /(NN-1);
	    fExpectedMBRVrmsCut->FixParameter(0, slope);
	    double exp_VrmsCut = fExpectedMBRVrmsCut->Eval(NeA/1e9);
	    double p_exp_cut = pow(exp_VrmsCut,2)/50. * 1e6;
	    double I_limit_cut = thPowerCut[ipeaks]->GetMean()/p_exp_cut;
	    tgMBRFluxLimitCut[ipeaks]->SetPoint(i,slope,I_limit_cut);


	  }

	if(ipeaks == 0) {
	  tgMBRFluxLimitCut[ipeaks]->SetLineColor(1);
	  tgMBRFluxLimitCut[ipeaks]->SetFillColor(1);
	
	} else {

	  tgMBRFluxLimitCut[ipeaks]->SetLineColor(ipeaks * 2);
	  tgMBRFluxLimitCut[ipeaks]->SetFillColor(ipeaks * 2);

	}


	tgMBRFluxLimitCut[ipeaks]->SetLineWidth(2);
	tgMBRFluxLimitCut[ipeaks]->SetTitle(Form("MBR Intensity at Ne/1e9 = %f",NeA/1e9));
	tgMBRFluxLimitCut[ipeaks]->GetYaxis()->SetTitle("I/I_{0} [I_{0}=4#times 10^{-16} W/m^{2}/Hz]");
	tgMBRFluxLimitCut[ipeaks]->GetXaxis()->SetTitle("alpha");

	tgMBRFluxLimitCut[ipeaks]->Write(Form("tgMBRFluxLimitCut_sgn%d_peak%d",ich,ipeaks));
	thPowerCut[ipeaks]->Write(Form("thPowerCut_sgn%d_peak%d",ich,ipeaks));
	

       }

  fTreeOut->Write();
 
  if(plot == 100) 
    {

      hh1_exp->Write();
      hh2_exp->Write();
      h_AeffCal->Write();
      h_Aeff->Write();
      h_Cal->Write();
    }
  

  tgCalib->Write("tgCalib");
  tgAntenna->Write("tgAntenna");
  tgAeff->Write("tgAeff");

  fout->Close();
}
 
void Put_Pointer_into_Vector(vector<double> *Vp, vector<double> & V){
  V.clear();
  for (unsigned int j=0; j<Vp->size(); ++j) 
    V.push_back(Vp->at(j));
}

void StripComments(string& line) {
  
  TString buffer(line);
  TPRegexp comments("#.*$");
  comments.Substitute(buffer, "");
  line = string(buffer.Data());
  
}

void AlignSignals(double Tc1, double Tc2, vector<double> *fVA, vector<double> *fTA, vector<double> *fVB, vector<double> *fTB, double fNeA, double fNeB, double & fdt_alignment, double & fchi2){
  //
  // Align in time fVB:fTB with respect to fVA:fTA
  // all the times are in ns
  // 
  int dt_bin = 0;
  fdt_alignment = 0.;
  fchi2 = 1e9;
  //
  int Nbins_t;
  double tmin,tmax,stept;
  Nbins_t = fTA->size(); tmin = fTA->at(0)/ns; tmax = fTA->at(Nbins_t-1)/ns;
  stept = fTA->at(1)/ns-fTA->at(0)/ns;

  // time range to find the alignment 
  //double Tc1 = -306;
  //double Tc2 = -296;
  Tc1 = Tc1/ns;
  Tc2 = Tc2/ns;

  //double Tc = (Tc1+Tc2)/2.;
  //int iTc = (Tc-tmin)/stept;
  //int diTc = (Tc-Tc1)/stept; 
  int iTc1 = (Tc1-tmin)/stept;
  int iTc2 = (Tc2-tmin)/stept;

  // how many ns do you want to shift the time trace to find the alignment 
  // note that 40GS/s --> 25 ps = 0.025 ns  ( 1 ns <-> 40 time bins)
  //double dtmax = 1.; //ns
  double dtmax = 2.; //ns
  int idtmax = dtmax/stept;

  //printf("step=%f: t =  %f - %f  <->  it = %d - %d   <->  shift_max %f (%d) \n",stept,Tc1,Tc2,iTc1,iTc2,dtmax,idtmax);  
  //mprintf("fNeA=%le   fNeB=%le \n",fNeA,fNeB);

  double chi2min = 1e9;
  for (int idt=-idtmax; idt<=idtmax; ++idt) {
    double chi2 = 0.;
    //printf("\nma");
    for (int it=iTc1; it<iTc2; ++it) {
      double y1 = fVA->at(it);
      double y2 = fVB->at(it+idt);
      y2 = y2 * fNeA/fNeB;
      //printf("%d  %f:  (%d,%d)   v0=%f  v1=%f \n",idt,fTA->at(it)/ns,iTc1,iTc2,y1*1e3,y2*1e3); 
      chi2 = chi2 + pow((y1-y2),2)/2./1e-3;
    }
    //if (EvtId==13813)
      //printf("t = %f - %f    idt=%d    chi2 = %f (%f) \n",Tc1,Tc2,idt,chi2,chi2min);
    if (chi2<chi2min) {
      chi2min = chi2;
      dt_bin = idt;
    }
  }

  //cout << fVA->at(iTc1) << " " << fVB->at(iTc1) << endl;

  fdt_alignment = dt_bin * stept * ns;
  fchi2 = chi2min;

  //cout << " " << fdt_alignment/ns << endl;

  //printf("\n");
  //printf("dt_bin=%d   dt=%f   chi2=%f \n",dt_bin,fdt_alignment/ns,fchi2);
  //printf("\n");

}



void TransformFFT(vector<double> *Vall, vector<double> *Tall, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq, vector<double> & F, vector<double> & P, vector<double> & VCut, vector<double> & T, vector<double> & PCut){
 
  double *in;
  fftw_complex *out_freq;
  fftw_plan p;	
  double *out_trace;
  fftw_plan p_trace;
  vector<double> Pmean;
  vector<double> V;
  T.clear();


  double ttmp=0;
  //select the trace between tmin and tmax
  for (unsigned int j=0; j<Tall->size(); ++j) {
    ttmp = Tall->at(j) - deltaT;
    if (ttmp<tmin || ttmp>tmax) continue;
    T.push_back(ttmp);
    V.push_back(Vall->at(j));
  }


  int N = V.size();
  Pmean.clear();
    
  in = (double*) fftw_malloc(sizeof(double) * N);
  out_freq = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N/2+1));
  p = fftw_plan_dft_r2c_1d(N, in, out_freq, FFTW_ESTIMATE);
    
  out_trace = (double*) fftw_malloc(sizeof(double) * N);
  p_trace = fftw_plan_dft_c2r_1d(N,out_freq, out_trace, FFTW_ESTIMATE);
  
  for(int i=0;i<N;i++) 
    in[i] = double(V[i]);
  
  fftw_execute(p);
    
  double dTime = fabs(T[1]-T[0]);
  double dFreq = 1./dTime/N;
  for(int i=0;i<(N/2+1);i++) {

    double spectrum = out_freq[i][0]*out_freq[i][0] + out_freq[i][1]*out_freq[i][1];

    F.push_back(i*dFreq);
    P.push_back(spectrum/N);
    bool cut=false;
    
    for (unsigned int f=0; f<MinFreq.size();f++)
      if(F[i]/1e9>MinFreq.at(f) && F[i]/1e9<MaxFreq.at(f)) 
	cut = true;

    if(!cut) 
      PCut.push_back(spectrum/N);
    else 
      PCut.push_back(1e-20);
	
    if(cut) {
      out_freq[i][0] = 0.;
      out_freq[i][1] = 0.;
    }
  }

  fftw_execute(p_trace);
  for(int i=0;i<N;i++) 
    VCut.push_back(out_trace[i]/N);

}

void TransformFFT2(vector<double> *Vall, vector<double> *Tall, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq, vector<double> & ReP, vector<double> & ImP)
{
  
  double *in;
  fftw_complex *out_freq;
  fftw_plan p;	
  double *out_trace;
  fftw_plan p_trace;
  vector<double> Pmean;
  vector<double> V, T;
  
  double ttmp=0;
  //select the trace between tmin and tmax
  for (unsigned int j=0; j<Tall->size(); ++j) {
    ttmp = Tall->at(j) - deltaT;
    if (ttmp<tmin || ttmp>tmax) continue;
    T.push_back(ttmp);
    V.push_back(Vall->at(j));
  }


  int N = V.size();
  Pmean.clear();
    
  //double dTime = fabs(T[1]-T[0]);
  //double dFreq = 1./dTime/N;

  in = (double*) fftw_malloc(sizeof(double) * N);
  out_freq = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N/2+1));
  p = fftw_plan_dft_r2c_1d(N, in, out_freq, FFTW_ESTIMATE);
    
  out_trace = (double*) fftw_malloc(sizeof(double) * N);
  p_trace = fftw_plan_dft_c2r_1d(N,out_freq, out_trace, FFTW_ESTIMATE);
  
  for(int i=0;i<N;i++) {
    in[i] = double(V[i]);
  }
  fftw_execute(p);
    
  for(int i=0;i<(N/2+1);i++) {
    //double spectrum = out_freq[i][0]*out_freq[i][0]+out_freq[i][1]*out_freq[i][1];

    ReP.push_back(out_freq[i][0]);
    ImP.push_back(out_freq[i][1]);
    
  }
}

void GetNe(vector<double> *sgn, vector<double> *time, double tped1, double tped2, double t1, double t2, double bunchlength, double &ne, double &ped, double &sgn_max, double &current){ 

  ne = 0.;
  ped = 0.;

  sgn_max = -1e9;
   
  //double ne = 0.; 
  //double ped = 0.;
  int jj_ped = 0; 
  int jj_sgn = 0; 

  int n = sgn->size();
  for (int j=0; j<n; j++) {
    double t = time->at(j);
    if (t>tped1 && t<tped2){
      ped = ped - sgn->at(j);
      jj_ped++;
    }
    if (t>t1 && t<t2){
      ne = ne - sgn->at(j);
      jj_sgn++;
      if (sgn_max < fabs(sgn->at(j)) ){
	sgn_max = fabs(sgn->at(j));
      }
    }
  }
  ped = ped / double (jj_ped);
  ne = ne - ped * double (jj_sgn);

  // Carica ottenuta dall'oscilloscopio
  // Q [C] = g Sgn[V*s] / Resistenza 
  // Numero di elettroni
  // Ne = Q[C] / qe
  double qe = 1.602176565 * 1e-19; // C
  double g = 10;  // calibrazione monitor bergoz Q(beam)/Q(out) @ 50Ohm = g = 10(btf)  
  double Resistenza = 50.; // Ohm
  double Sgn = ne * ( time->at(11) - time->at(10) ); // questo e' Sum_i Segnale_i dt [V sec]
  ne = g * Sgn / Resistenza / qe;

  // Corrente per bunch in mA
  // I [mA] = g Sgn[pV*s] / Resistenza / bunchlength [ns] 
  current = g * Sgn*1e12 / Resistenza / bunchlength;
  //printf("%le   %f \n ",ne,current);

  //printf("t1=%f   t2=%f    Ne=%le \n",t1/ns,t2/ns,ne);

}



double GetVrms_bin(vector<double> *sgn, vector<double> *time, double deltaT, double t1, double t2){ 

  // nuovo metodo per ottenere il VRMS 
  // e' piu' rapido in quanto si fa il loop solo sui bin definiti entro t1 e t2
  // e' da usare per fare delle moving average
  // E' stato testato e' con i1 e i2 definiti in questo modo da lo stesso 
  // VRMS di GetVrms

  // RMS
  int Nav = 0;
  double Vrms = 0.; 
  int N = sgn->size();

  double tmin = time->at(0)/ns;
  double step = time->at(1)/ns - time->at(0)/ns;
  // nota1:  in GetVrms si ha t-deltat > t1  e quindi qui bisogna avere t > t1 + deltat  
  int i1 = ( t1/ns+deltaT/ns - tmin )/step + 1; 
  int i2 = ( t2/ns+deltaT/ns - tmin )/step + 1; 

  //printf("i1=%d  t1=%f   i2=%d  t2=%f \n",i1,t1/ns,i2,t2/ns);

  if (i1<0) i1 = 0;
  if (i2<0) i2 = 0;
  if (i1>N) i1 = N;
  if (i2>N) i2 = N;
  if (i1>i2) i1 = i2;
  //if (i2<=i1) printf("i1=%d  t1=%f   i2=%d  t2=%f \n",i1,t1/ns,i2,t2/ns);

  //if (DEBUG) printf("i1=%d  i2=%d \n",i1,i2);

  for (int i=i1; i<i2; i++) { 
    Vrms += pow(sgn->at(i),2);
    Nav++;
  }
  if (Nav>0) 
    Vrms = sqrt(Vrms/double(Nav));

  /**
   // Calcolo in GetVrms 

  for (int i=0; i<N; i++) { 
    double t = time->at(i)-deltaT;
    if (t>t1 && t<t2){ 
      Vrms += pow(sgn->at(i),2);
      Nav++;
    }
  }
  Vrms = sqrt(Vrms/double(Nav));
  */

  return Vrms;
}

double GetVrms(vector<double> *sgn, vector<double> *time, double deltaT, double t1, double t2){ 

  // RMS
  int Nav = 0;
  double Vrms = 0.; 
  int N = sgn->size();

  for (int i=0; i<N; i++) { 
    double t = time->at(i)-deltaT;
    if (t>t1 && t<t2){ 
      //if (DEBUG) printf("i=%d \n",i);
      Vrms += pow(sgn->at(i),2);
      Nav++;
    }
  }
  Vrms = sqrt(Vrms/double(Nav));

  return Vrms;
}

double GetVrms1(vector<double> sgn, vector<double> time, double deltaT, double t1, double t2){ 

  // RMS
  int Nav = 0;
  double Vrms = 0.; 
  int N = sgn.size();
  for (int i=0; i<N; i++) { 
    double t = time.at(i)-deltaT;
    if (t>t1 && t<t2){ 
      Vrms += pow(sgn.at(i),2);
      Nav++;
    }
  }
  Vrms = sqrt(Vrms/double(Nav));

  return Vrms;
}


double PowerFromSpectrum(vector<double> *P, vector<double> *F, vector<double> MinFreq, vector<double> MaxFreq){
  //
  // return the power in W
  //
  double yy = 0.;
  for(int i=0;i<F->size();i++) {
    
    bool cut = false;
    for (unsigned int f=0; f<MinFreq.size();f++){
      if ( F->at(i)/1e9 > MinFreq.at(f) && 
           F->at(i)/1e9 < MaxFreq.at(f) )
	cut = true; 
    }

    if (!cut)
      yy = yy + P->at(i);

  }

  return ( yy / R / F->size() );

}

double PowerFromSpectrum1(vector<double> P, vector<double> F, vector<double> MinFreq, vector<double> MaxFreq){
  //
  // return the power in W
  //
  double yy = 0.;
  for(int i=0;i<F.size();i++) {
    
    bool cut = false;
    for (unsigned int f=0; f<MinFreq.size();f++){
      if ( F.at(i)/1e9 > MinFreq.at(f) && 
           F.at(i)/1e9 < MaxFreq.at(f) )
	cut = true; 
    }

    if (!cut)
      yy = yy + P.at(i);

  }

  return ( yy / R / F.size() );

}


void BookDSTData(){

  fData.SetBranchAddress("RunId",&RunId);
  fData.SetBranchAddress("EvtId",&EvtId);
  fData.SetBranchAddress("T1",&T1);
  fData.SetBranchAddress("T2",&T2);
  fData.SetBranchAddress("T3",&T3);
  fData.SetBranchAddress("T4",&T4);
  fData.SetBranchAddress("Ch1",&Ch1);
  fData.SetBranchAddress("Ch2",&Ch2);
  fData.SetBranchAddress("Ch3",&Ch3);
  fData.SetBranchAddress("Ch4",&Ch4);

  if(IsDecember2012)
    fData.SetBranchAddress("FixVertGain",&FixVertGain[0]);


}

void InitOutputTree(){

  evtid = 0;
  runid = 0.;
  bunch_length = 0.;
  target = 0;
  Ne = 0.; Ne1 = 0.; Ne2 = 0.; Ne3 = 0.; Ne4 = 0.;
  Current = 0.;
  ped = 0.; t_beam_max = 0.; Ch1_max = 0.;

  chi2_al = 0.;
  dt_al = 0.;

  for (int ipeak=0; ipeak<NPeak; ipeak++) { 
    vrms_sgn_ch1[ipeak] = 0.; vrms_sgn_ch2[ipeak] = 0.; vrms_sgn_ch3[ipeak] = 0.; vrms_sgn_ch4[ipeak] = 0.; 
    vrms_bkg_ch1[ipeak] = 0.; vrms_bkg_ch2[ipeak] = 0.; vrms_bkg_ch3[ipeak] = 0.; vrms_bkg_ch4[ipeak] = 0.; 
    vrms_sgn_ch1_fsel1[ipeak] = 0.; vrms_sgn_ch2_fsel1[ipeak] = 0.; vrms_sgn_ch3_fsel1[ipeak] = 0.; vrms_sgn_ch4_fsel1[ipeak] = 0.; 
    vrms_bkg_ch1_fsel1[ipeak] = 0.; vrms_bkg_ch2_fsel1[ipeak] = 0.; vrms_bkg_ch3_fsel1[ipeak] = 0.; vrms_bkg_ch4_fsel1[ipeak] = 0.; 
    vrms_sgn_ch1_fsel2[ipeak] = 0.; vrms_sgn_ch2_fsel2[ipeak] = 0.; vrms_sgn_ch3_fsel2[ipeak] = 0.; vrms_sgn_ch4_fsel2[ipeak] = 0.; 
    vrms_bkg_ch1_fsel2[ipeak] = 0.; vrms_bkg_ch2_fsel2[ipeak] = 0.; vrms_bkg_ch3_fsel2[ipeak] = 0.; vrms_bkg_ch4_fsel2[ipeak] = 0.; 

    /*
    vrms_exp1_sgn_ch1[ipeak] = 0.; vrms_exp1_sgn_ch2[ipeak] = 0.; vrms_exp1_sgn_ch3[ipeak] = 0.; vrms_exp1_sgn_ch4[ipeak] = 0.; 
    vrms_exp1_sgn_ch1_fsel1[ipeak] = 0.; vrms_exp1_sgn_ch2_fsel1[ipeak] = 0.; vrms_exp1_sgn_ch3_fsel1[ipeak] = 0.; vrms_exp1_sgn_ch4_fsel1[ipeak] = 0.; 
    vrms_exp1_sgn_ch1_fsel2[ipeak] = 0.; vrms_exp1_sgn_ch2_fsel2[ipeak] = 0.; vrms_exp1_sgn_ch3_fsel2[ipeak] = 0.; vrms_exp1_sgn_ch4_fsel2[ipeak] = 0.; 
    
    vrms_exp2_sgn_ch1[ipeak] = 0.; vrms_exp2_sgn_ch2[ipeak] = 0.; vrms_exp2_sgn_ch3[ipeak] = 0.; vrms_exp2_sgn_ch4[ipeak] = 0.; 
    vrms_exp2_sgn_ch1_fsel1[ipeak] = 0.; vrms_exp2_sgn_ch2_fsel1[ipeak] = 0.; vrms_exp2_sgn_ch3_fsel1[ipeak] = 0.; vrms_exp2_sgn_ch4_fsel1[ipeak] = 0.; 
    vrms_exp2_sgn_ch1_fsel2[ipeak] = 0.; vrms_exp2_sgn_ch2_fsel2[ipeak] = 0.; vrms_exp2_sgn_ch3_fsel2[ipeak] = 0.; vrms_exp2_sgn_ch4_fsel2[ipeak] = 0.; 
    */

    vrms_dedx_exp1_sgn_ch1[ipeak] = 0.; vrms_dedx_exp1_sgn_ch2[ipeak] = 0.; vrms_dedx_exp1_sgn_ch3[ipeak] = 0.; vrms_dedx_exp1_sgn_ch4[ipeak] = 0.; 
    vrms_dedx_exp1_sgn_ch1_fsel1[ipeak] = 0.; vrms_dedx_exp1_sgn_ch2_fsel1[ipeak] = 0.; vrms_dedx_exp1_sgn_ch3_fsel1[ipeak] = 0.; vrms_dedx_exp1_sgn_ch4_fsel1[ipeak] = 0.; 
    vrms_dedx_exp1_sgn_ch1_fsel2[ipeak] = 0.; vrms_dedx_exp1_sgn_ch2_fsel2[ipeak] = 0.; vrms_dedx_exp1_sgn_ch3_fsel2[ipeak] = 0.; vrms_dedx_exp1_sgn_ch4_fsel2[ipeak] = 0.; 
    
    vrms_dedx_exp2_sgn_ch1[ipeak] = 0.; vrms_dedx_exp2_sgn_ch2[ipeak] = 0.; vrms_dedx_exp2_sgn_ch3[ipeak] = 0.; vrms_dedx_exp2_sgn_ch4[ipeak] = 0.; 
    vrms_dedx_exp2_sgn_ch1_fsel1[ipeak] = 0.; vrms_dedx_exp2_sgn_ch2_fsel1[ipeak] = 0.; vrms_dedx_exp2_sgn_ch3_fsel1[ipeak] = 0.; vrms_dedx_exp2_sgn_ch4_fsel1[ipeak] = 0.; 
    vrms_dedx_exp2_sgn_ch1_fsel2[ipeak] = 0.; vrms_dedx_exp2_sgn_ch2_fsel2[ipeak] = 0.; vrms_dedx_exp2_sgn_ch3_fsel2[ipeak] = 0.; vrms_dedx_exp2_sgn_ch4_fsel2[ipeak] = 0.; 
    
  }   
  
}


void BookOutputTree(TTree *fTreeOut){

  fTreeOut->Branch("runid",&runid,"runid/I");
  fTreeOut->Branch("evtid",&evtid,"evtid/I");
  fTreeOut->Branch("Ne",&Ne,"Ne/D");
  fTreeOut->Branch("Ne1",&Ne1,"Ne1/D");
  fTreeOut->Branch("Ne2",&Ne2,"Ne2/D");
  fTreeOut->Branch("Ne3",&Ne3,"Ne3/D");
  fTreeOut->Branch("Ne4",&Ne4,"Ne4/D");
  fTreeOut->Branch("ped",&ped,"ped/D");
  fTreeOut->Branch("Current",&Current,"Current/D");
  fTreeOut->Branch("chi2_al",&chi2_al,"chi2_al/D");
  fTreeOut->Branch("dt_al",&dt_al,"dt_al/D");
 
 // WARNING: il [3] deve essere [NPeak]: usa sprintf
  fTreeOut->Branch("vrms_sgn_ch1",&vrms_sgn_ch1,"vrms_sgn_ch1[3]/D");
  fTreeOut->Branch("vrms_sgn_ch2",&vrms_sgn_ch2,"vrms_sgn_ch2[3]/D");
  fTreeOut->Branch("vrms_sgn_ch3",&vrms_sgn_ch3,"vrms_sgn_ch3[3]/D");
  fTreeOut->Branch("vrms_sgn_ch4",&vrms_sgn_ch4,"vrms_sgn_ch4[3]/D");

  fTreeOut->Branch("vrms_bkg_ch1",&vrms_bkg_ch1,"vrms_bkg_ch1[3]/D");
  fTreeOut->Branch("vrms_bkg_ch2",&vrms_bkg_ch2,"vrms_bkg_ch2[3]/D");
  fTreeOut->Branch("vrms_bkg_ch3",&vrms_bkg_ch3,"vrms_bkg_ch3[3]/D");
  fTreeOut->Branch("vrms_bkg_ch4",&vrms_bkg_ch4,"vrms_bkg_ch4[3]/D");

  fTreeOut->Branch("vrms_sgn_ch1_fsel1",&vrms_sgn_ch1_fsel1,"vrms_sgn_ch1_fsel1[3]/D");
  fTreeOut->Branch("vrms_sgn_ch2_fsel1",&vrms_sgn_ch2_fsel1,"vrms_sgn_ch2_fsel1[3]/D");
  fTreeOut->Branch("vrms_sgn_ch3_fsel1",&vrms_sgn_ch3_fsel1,"vrms_sgn_ch3_fsel1[3]/D");
  fTreeOut->Branch("vrms_sgn_ch4_fsel1",&vrms_sgn_ch4_fsel1,"vrms_sgn_ch4_fsel1[3]/D");
  fTreeOut->Branch("vrms_bkg_ch1_fsel1",&vrms_bkg_ch1_fsel1,"vrms_bkg_ch1_fsel1[3]/D");
  fTreeOut->Branch("vrms_bkg_ch2_fsel1",&vrms_bkg_ch2_fsel1,"vrms_bkg_ch2_fsel1[3]/D");
  fTreeOut->Branch("vrms_bkg_ch3_fsel1",&vrms_bkg_ch3_fsel1,"vrms_bkg_ch3_fsel1[3]/D");
  fTreeOut->Branch("vrms_bkg_ch4_fsel1",&vrms_bkg_ch4_fsel1,"vrms_bkg_ch4_fsel1[3]/D");

  /*
  fTreeOut->Branch("vrms_exp1_sgn_ch1",&vrms_exp1_sgn_ch1,"vrms_exp1_sgn_ch1[3]/D");
  fTreeOut->Branch("vrms_exp1_sgn_ch2",&vrms_exp1_sgn_ch2,"vrms_exp1_sgn_ch2[3]/D");
  fTreeOut->Branch("vrms_exp1_sgn_ch3",&vrms_exp1_sgn_ch3,"vrms_exp1_sgn_ch3[3]/D");
  fTreeOut->Branch("vrms_exp1_sgn_ch4",&vrms_exp1_sgn_ch4,"vrms_exp1_sgn_ch4[3]/D");

  fTreeOut->Branch("vrms_exp1_sgn_ch1_fsel1",&vrms_exp1_sgn_ch1_fsel1,"vrms_exp1_sgn_ch1_fsel1[3]/D");
  fTreeOut->Branch("vrms_exp1_sgn_ch2_fsel1",&vrms_exp1_sgn_ch2_fsel1,"vrms_exp1_sgn_ch2_fsel1[3]/D");
  fTreeOut->Branch("vrms_exp1_sgn_ch3_fsel1",&vrms_exp1_sgn_ch3_fsel1,"vrms_exp1_sgn_ch3_fsel1[3]/D");
  fTreeOut->Branch("vrms_exp1_sgn_ch4_fsel1",&vrms_exp1_sgn_ch4_fsel1,"vrms_exp1_sgn_ch4_fsel1[3]/D");
 
  fTreeOut->Branch("vrms_exp2_sgn_ch1",&vrms_exp2_sgn_ch1,"vrms_exp2_sgn_ch1[3]/D");
  fTreeOut->Branch("vrms_exp2_sgn_ch2",&vrms_exp2_sgn_ch2,"vrms_exp2_sgn_ch2[3]/D");
  fTreeOut->Branch("vrms_exp2_sgn_ch3",&vrms_exp2_sgn_ch3,"vrms_exp2_sgn_ch3[3]/D");
  fTreeOut->Branch("vrms_exp2_sgn_ch4",&vrms_exp2_sgn_ch4,"vrms_exp2_sgn_ch4[3]/D");
 
  fTreeOut->Branch("vrms_exp2_sgn_ch1_fsel1",&vrms_exp2_sgn_ch1_fsel1,"vrms_exp2_sgn_ch1_fsel1[3]/D");
  fTreeOut->Branch("vrms_exp2_sgn_ch2_fsel1",&vrms_exp2_sgn_ch2_fsel1,"vrms_exp2_sgn_ch2_fsel1[3]/D");
  fTreeOut->Branch("vrms_exp2_sgn_ch3_fsel1",&vrms_exp2_sgn_ch3_fsel1,"vrms_exp2_sgn_ch3_fsel1[3]/D");
  fTreeOut->Branch("vrms_exp2_sgn_ch4_fsel1",&vrms_exp2_sgn_ch4_fsel1,"vrms_exp2_sgn_ch4_fsel1[3]/D");
  */

  fTreeOut->Branch("vrms_dedx_exp1_sgn_ch1",&vrms_dedx_exp1_sgn_ch1,"vrms_dedx_exp1_sgn_ch1[3]/D");
  fTreeOut->Branch("vrms_dedx_exp1_sgn_ch2",&vrms_dedx_exp1_sgn_ch2,"vrms_dedx_exp1_sgn_ch2[3]/D");
  fTreeOut->Branch("vrms_dedx_exp1_sgn_ch3",&vrms_dedx_exp1_sgn_ch3,"vrms_dedx_exp1_sgn_ch3[3]/D");
  fTreeOut->Branch("vrms_dedx_exp1_sgn_ch4",&vrms_dedx_exp1_sgn_ch4,"vrms_dedx_exp1_sgn_ch4[3]/D");

  fTreeOut->Branch("vrms_dedx_exp1_sgn_ch1_fsel1",&vrms_dedx_exp1_sgn_ch1_fsel1,"vrms_dedx_exp1_sgn_ch1_fsel1[3]/D");
  fTreeOut->Branch("vrms_dedx_exp1_sgn_ch2_fsel1",&vrms_dedx_exp1_sgn_ch2_fsel1,"vrms_dedx_exp1_sgn_ch2_fsel1[3]/D");
  fTreeOut->Branch("vrms_dedx_exp1_sgn_ch3_fsel1",&vrms_dedx_exp1_sgn_ch3_fsel1,"vrms_dedx_exp1_sgn_ch3_fsel1[3]/D");
  fTreeOut->Branch("vrms_dedx_exp1_sgn_ch4_fsel1",&vrms_dedx_exp1_sgn_ch4_fsel1,"vrms_dedx_exp1_sgn_ch4_fsel1[3]/D");
 
  fTreeOut->Branch("vrms_dedx_exp2_sgn_ch1",&vrms_dedx_exp2_sgn_ch1,"vrms_dedx_exp2_sgn_ch1[3]/D");
  fTreeOut->Branch("vrms_dedx_exp2_sgn_ch2",&vrms_dedx_exp2_sgn_ch2,"vrms_dedx_exp2_sgn_ch2[3]/D");
  fTreeOut->Branch("vrms_dedx_exp2_sgn_ch3",&vrms_dedx_exp2_sgn_ch3,"vrms_dedx_exp2_sgn_ch3[3]/D");
  fTreeOut->Branch("vrms_dedx_exp2_sgn_ch4",&vrms_dedx_exp2_sgn_ch4,"vrms_dedx_exp2_sgn_ch4[3]/D");
  
  fTreeOut->Branch("vrms_dedx_exp2_sgn_ch1_fsel1",&vrms_dedx_exp2_sgn_ch1_fsel1,"vrms_dedx_exp2_sgn_ch1_fsel1[3]/D");
  fTreeOut->Branch("vrms_dedx_exp2_sgn_ch2_fsel1",&vrms_dedx_exp2_sgn_ch2_fsel1,"vrms_dedx_exp2_sgn_ch2_fsel1[3]/D");
  fTreeOut->Branch("vrms_dedx_exp2_sgn_ch3_fsel1",&vrms_dedx_exp2_sgn_ch3_fsel1,"vrms_dedx_exp2_sgn_ch3_fsel1[3]/D");
  fTreeOut->Branch("vrms_dedx_exp2_sgn_ch4_fsel1",&vrms_dedx_exp2_sgn_ch4_fsel1,"vrms_dedx_exp2_sgn_ch4_fsel1[3]/D");
 


}


bool SelTriggers(int fRunId, int fEvtId, double fFixVertGain[4], 
		 vector<double> *ch1,  vector<double> *t1, 
		 vector<double> *ch2,  vector<double> *t2, 
		 vector<double> *ch3,  vector<double> *t3, 
		 vector<double> *ch4,  vector<double> *t4,
		 double ne_min, double ne_max, double tb1, double tb2, double ts1, double ts2,
		 bool isdecember2012,
		 bool selch[5]
		 ) {

  bool fSel = 1;

  if ( t1->size()>1 && t1->size()==t2->size() && t1->size()==t3->size() && t1->size()==t4->size() ){
    // 
  } else {
    fSel = 0;
  }

  // Intenisty_Scan_3ns_C2_CO_T0
  if (fRunId==12132015 || fRunId==12131934 || fRunId==12132031)
    fSel = 0;


  // 1ns-T4-Cross-Dec2012
  if (fRunId==12161154 && fEvtId<=42169)
    fSel = 0;

  // 1ns-T5-Cross-Dec2012
  if (fRunId==12161115 && fEvtId>41662)
    fSel = 0;



  // Intenisty_Scan_1ns_C2_cross_T0
  /**
  if (fEvtId>77150 && fEvtId<77400) fSel = 0;
  if (fEvtId>77740 && fEvtId<77840) fSel = 0;
  if (fEvtId>78130 && fEvtId<78400) fSel = 0;
  if (fEvtId>78370 && fEvtId<78690) fSel = 0;
  */

  // Intenisty_Scan_1ns_C2_cross_T3
  // in questi trigger e' cambiata la lunghezza del bunch (da 1.5 ns a 3 ns)
  // if (fRunId==12151302 && fEvtId>=17756 && fEvtId<=17783) fSel = 0;

  // Intenisty_Scan_1ns_C2_copo_T2
  if (fRunId==12150840 && fEvtId>12260 && fEvtId<12310) fSel = 0;
  if (fRunId==12150840 && fEvtId>12000 && fEvtId<12040) fSel = 0;

// Intenisty_Scan_1ns_C2_copo_T3
  if (fRunId==12150828 && fEvtId>11500) fSel = 0;

// Intenisty_Scan_1ns_C2_copo_T4
  if (fRunId==12150714 && fEvtId>10940 && fEvtId<10970) fSel = 0;

// Intenisty_Scan_3ns_C2_cross_T3
  // if (fRunId==12110539 && fEvtId<23050) fSel = 0;


  if (fSel==1){

    double BunchLength_dummy = 1.;
    double fne, fped, fch1_max, fCurrent;
    GetNe(ch1, t1, tb1, tb2, ts1, ts2, BunchLength_dummy, fne, fped, fch1_max, fCurrent);
    if (fne<ne_min || fne>ne_max) fSel = 0;
    
    if (isdecember2012){
      if (IsSaturated(fRunId,fEvtId,1,ch1,fFixVertGain[0])) fSel = 0;
      if (selch[2] && IsSaturated(fRunId,fEvtId,2,ch2,fFixVertGain[2-1])) fSel = 0;
      if (selch[3] && IsSaturated(fRunId,fEvtId,3,ch3,fFixVertGain[3-1])) fSel = 0;
      if (selch[4] && IsSaturated(fRunId,fEvtId,4,ch4,fFixVertGain[4-1])) fSel = 0;
    }

  }

  return (fSel);
}


bool IsSaturated(int fRunId, int fEvtId, int ich, vector<double> *ch, double fixvertgain){
  bool sat = false;
  int nbins_sat = 0;
  int n = ch->size();
  for (int j=0; j<n; j++) {
    if (TMath::Abs(ch->at(j)/fixvertgain) >= 4){
      sat = true;
      nbins_sat++;
    }
  }

  //sat = false;

  if (sat)
    cout << "WARNING: channel "<< ich << " has " << nbins_sat << " bins saturated. Skip this trigger: " << fRunId << " " << fEvtId << endl;
  return (sat);
}

double Aeff(double G, double nu){
  // return antenna effective area in m2
  double A = pow(10.,G/10.)/(4.*TMath::Pi()) * pow(SpeedLight/(nu*GHz),2) ;
  return (A);
}

double Calib_Dec2012(double nu){
  // return horn gain in dBi
  double y = tgCalib->Eval(nu*1e9,0,"");
  return (pow(10,y/10.)); 
}

double ExpectedVrms(double ne, double alpha, int antenna, double BunchLength, int selfreq,vector<double> *time, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq){
  
  vector<double> T;
  
  double ttmp=0;
  //select the trace between tmin and tmax
  for (unsigned int j=0; j<time->size(); ++j) {
    ttmp = time->at(j) - deltaT;
    if (ttmp<tmin || ttmp>tmax) continue;
    T.push_back(ttmp);
  }

 
  int N = T.size();
  // cout<<" Tsize"<< N<<endl;
  // flux density 

  // Ghoram
  double I0 = 4.*1e-16; // W/m2/Hz
  //I0 = I0 * 1e-4;
  double E0 = 3.36 *1e17; // eV

  // test beam
  //double ne = 1e9;
  double E = ne * 510.*1e6; // eV

  // factor to account the different geometries in the two experiments
  double g = 0.1*(10./BunchLength);
  // se il bunch e' 3 ns moltiplicare per 10./3
  // se il bunch e' 1,5 ns moltiplicare per 10./1,5

  // expected density flux
  double I = I0 * g * pow(E/E0,alpha);

  //if(ne/1e9 > 0.9 && ne/1e9 < 1.1 && alpha == 2 && selfreq == 0)
  //if(alpha == 2 && selfreq == 0)
  //printf("Expected VRMS: %f %f \n", ne/1e9, g * pow(E/E0,alpha));

  // expected signal 
  double P     = 0;
  double P_cal = 0;
  double P_dbm = 0;
  double P_cal_dbm = 0;
  double V = 0;
  double V_cal = 0;

  
  double dTime = fabs(T[1]-T[0]);
  double dFreq = 1./dTime/double(N);
  double StepNu=dFreq;

  //for (int l=0;l<NNu;l++) {
  for (unsigned int l=1; l<(N/2+1);l++){
    double nu, lambda;
    //nu = NuMin + l*StepNu + StepNu/2.;

    nu = l*dFreq;
    //cout<<"l" <<l<<" "<<dFreq<<" " <<nu<<" " <<endl;

    lambda = SpeedLight/(nu*GHz); // m
    double gain=tgAntenna->Eval(nu);

    double dP=0;
    if(selfreq==1)
      {
	dP=0;
      }else{
      
      //if (antenna == 1) gain = G_LogPer(nu);
      //if (antenna == 2) gain = G_Horn(nu);

      dP = I * Aeff(gain,nu/GHz) * StepNu * Calib_Dec2012(nu/GHz);
     
    }
    P_cal = P_cal + dP;
    P     = P     + dP/Calib_Dec2012(nu/GHz);
    //cout<<"l" <<l<<" "<<dFreq/GHz<<" " <<nu/GHz<<" " <<gain<<" "<<Aeff(gain,nu/GHz)<<" "<<Calib_Dec2012(nu/GHz)<<" "<<StepNu/GHz<< " dP"<<dP<<" "<<P_cal<<endl;
  }
  // In order to conserve the total power, multiply all frequencies that occur in both sets — the positive and negative frequencies — by a factor of 2

  P_cal_dbm = 10.*log10(2*P_cal*1e3);
  V_cal = sqrt( 2*P_cal * R );

  P_dbm = 10.*log10(2*P*1e3);
  V = sqrt( 2*P * R );

  return (V_cal);

}


double ExpectedVrms1(double ne, double alpha, int antenna, double BunchLength, int selfreq,vector<double> time, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq){
  
  vector<double> T;
  
  double ttmp=0;
  //select the trace between tmin and tmax
  for (unsigned int j=0; j<time.size(); ++j) {
    ttmp = time.at(j) - deltaT;
    if (ttmp<tmin || ttmp>tmax) continue;
    T.push_back(ttmp);
  }


  int N = T.size();

  // flux density 

  // Ghoram
  double I0 = 4.*1e-16; // W/m2/Hz
  //I0 = I0 * 1e-4;
  double E0 = 3.36 *1e17; // eV

  // test beam
  //double ne = 1e9;
  double E = ne * 510.*1e6; // eV

  // factor to account the different geometries in the two experiments

  double g = 0.1 * (10./BunchLength);
  // se il bunch e' 3 ns moltiplicare per 10./3

  // expected density flux
  double I = I0 * g * pow(E/E0,alpha);

  //printf("%f \n", g * pow(E/E0,alpha));

  // expected signal 
  double P     = 0;
  double P_cal = 0;
  double P_dbm = 0;
  double P_cal_dbm = 0;
  double V = 0;
  double V_cal = 0;

  
  double dTime = fabs(T[1]-T[0]);
  double dFreq = 1./dTime/N;
  double StepNu=dFreq;



//for (int l=0;l<NNu;l++) {
  for (unsigned int l=1; l<(N/2+1);l++){
    double nu, lambda;
    //nu = NuMin + l*StepNu + StepNu/2.;

    nu = l*dFreq;
    //cout<<dTime<< " l  " <<l<<" "<<dFreq<<" " <<nu<<" " <<endl;

    lambda = SpeedLight/(nu*GHz); // m
    double gain=tgAntenna->Eval(nu);

    double dP=0;
    bool cut=false;
    
    for (unsigned int f=0; f<MinFreq.size();f++){
      //cout<<nu/GHz<<" "<<MinFreq.at(f)<<" "<<MaxFreq.at(f)<<endl;
      if(nu/GHz>MinFreq.at(f) && nu/GHz<MaxFreq.at(f)) 
	cut = true;
    }
    if(cut==true &&selfreq==1)
      {
	dP=0;
      }else{
      
      //if (antenna == 1) gain = G_LogPer(nu);
      //if (antenna == 2) gain = G_Horn(nu);

      dP = I * Aeff(gain,nu/GHz) * StepNu * Calib_Dec2012(nu/GHz);
     
    }
    P_cal = P_cal + dP;
    P     = P     + dP/Calib_Dec2012(nu/GHz);
    //cout<<"l " <<l<<" "<<dFreq/GHz<<" " <<nu/GHz<<" " <<gain<<" "<<Aeff(gain,nu/GHz)<<" "<<Calib_Dec2012(nu/GHz)<<" "<<StepNu/GHz<< " dP"<<dP<<" "<<P_cal<<endl;
  }
  // In order to conserve the total power, multiply all frequencies that occur in both sets — the positive and negative frequencies — by a factor of 2
  P_cal_dbm = 10.*log10(2*P_cal*1e3);
  V_cal = sqrt(2* P_cal * R );

  P_dbm = 10.*log10(P*1e3);
  V = sqrt(2* P * R );

  return (V_cal);

}

/*
  //for (int l=0;l<NNu;l++) {
  for (unsigned int l=0; l<N;l++){
    double nu, lambda;
    //nu = NuMin + l*StepNu + StepNu/2.;

    nu = l*dFreq; 
    lambda = SpeedLight/(nu*GHz); // m

    bool cut=false;
    
 for (unsigned int f=0; f<MinFreq.size();f++)
      if(nu*GHz>MinFreq.at(f) && nu*GHz/1e9<MaxFreq.at(f)) 
	cut = true;
    double dP=0;
    if(cut==true &&selfreq==1)
      {
	dP=0;
      }else{
      double gain=0;
      //if (antenna == 1) gain = G_LogPer(nu);
      //if (antenna == 2) gain = G_Horn(nu);
      gain=tgAntenna->Eval(nu*GHz);
      dP = I * Aeff(gain,nu) * StepNu*GHz * Calib_Dec2012(nu);
    }
    P_cal = P_cal + dP;
    P     = P     + dP/Calib_Dec2012(nu);
  }

  P_cal_dbm = 10.*log10(P_cal*1e3);
  V_cal = sqrt( P_cal * R );

  P_dbm = 10.*log10(P*1e3);
  V = sqrt( P * R );

  return (V_cal);

}
*/
  
//Gonzalo 10/11/2013
double ExpectedVrms_MBRdEdX(double ne, double alpha, int antenna, int Targets, double BunchLength, int selfreq,vector<double> *time, double deltaT, double tmin, double tmax, vector<double> MinFreq, vector<double> MaxFreq) {
  

vector<double> T;
  
  double ttmp=0;
  //select the trace between tmin and tmax
  for (unsigned int j=0; j<time->size(); ++j) {
    ttmp = time->at(j) - deltaT;
    if (ttmp<tmin || ttmp>tmax) continue;
    T.push_back(ttmp);
  }


  //int N = T.size();


  // Ghoram
  double I0 = 4.0*1e-16; // W/m2/Hz
  //I0 = I0 * 1e-4;
  //double E0 = 3.36 *1e17; // eV

  // test beam
  //double ne = 1e9;
  //double E = ne * 510.*1e6; // eV

  // factor to account the different geometries in the two experiments
  // double g = 1;//0.1 * (10. / BunchLenght);
  // se il bunch e' 3 ns moltiplicare per 10./3
  // se il bunch e' 1,5 ns moltiplicare per 10./1,5

  // expected density flux

  //double I = I0 * g * pow(E/E0,alpha);

  double Ng = 1.20e7; // Gorham number of electrons.
  //double ratio = FrascatidEdX[Target]/SLACdEdX[5];
  double ratio = FrascatidEdXFoV[Target]/SLACdEdXFoV[5];
  
  double Ls = 0.5;
  double Lf = 0.9;
  //double g = pow(Ls/Lf,2) * (1.3/3.0) * (10.0/BunchLength);
  double g = pow(Ls/Lf,2);
  double I = I0 * g * pow(ratio * ne/Ng , alpha);


  //if(alpha==2 && selfreq == 0)
  //cout << ne/1e9  << "\t"<< g * pow(ratio * ne/Ng , alpha) << endl;

  //printf("%f \n", g * pow(E/E0,alpha));

  // expected signal 
  double P     = 0;
  double P_cal = 0;
  double P_dbm = 0;
  double P_cal_dbm = 0;
  double V = 0;
  double V_cal = 0;

  
 // frequencies in GHz
  
  double bandwidth = 10. *1e-3; // GHz
  double NuMin = 2.0 - bandwidth/2.; //GHz
  double NuMax = 20.0 + bandwidth/2.; //GHz
  int NNu = (NuMax - NuMin)/bandwidth;  
  double StepNu = bandwidth * GHz;  //Hz
  

  //double dTime = fabs(T[1]-T[0]);
  //double dFreq = 1./dTime/N;
  //StepNu=dFreq; //Hz
  
 //bandwidth = StepNu / GHz; //GHz
  //NuMin = 2.0 - bandwidth/2.; //GHz
  //NuMax = 20.0 + bandwidth/2.; //GHz
  //NNu = (NuMax - NuMin)/bandwidth;  
  

  if(antenna == 3) {
    
    hh1_exp = new TH1D("hh1_exp","Power with calibration",NNu, NuMin, NuMax);
    hh2_exp = new TH1D("hh2_exp","Power without calibration",NNu, NuMin, NuMax);

    h_AeffCal = new TH1D("h_AeffCal","Aeff * Calibration",NNu, NuMin, NuMax);
    h_Aeff = new TH1D("h_Aeff","Aeff",NNu, NuMin, NuMax);
    h_Cal = new TH1D("h_Cal","Calibration",NNu, NuMin, NuMax);

  }

  for (int l=0;l<NNu;l++) {
    
    //for (unsigned int l=1; l<(N/2+1);l++){
    

    //double nu = (NuMin + l*StepNu + StepNu/2.) * GHz;

    double nu = (NuMin + l * StepNu/GHz) * GHz; // nu in Hz
    double gain = tgAntenna->Eval(nu); //nu in Hz

    //cout <<l << "\t"<< NuMin <<  "\t" << StepNu/GHz << "\t"  <<  nu << "\t gain: " << gain << "\t Aeff: "<< Aeff(gain,nu/GHz) << "\t calib: "<<  Calib_Dec2012(nu/GHz) << endl;
    
    bool cut=false;
    for (unsigned int f=0; f<MinFreq.size();f++){
      if(nu/GHz>MinFreq.at(f) && nu/GHz<MaxFreq.at(f))  //nu in GHz
	cut = true;
    }

    double dP=0;
    if(selfreq==1 && cut==true)
      {
	dP=0;
	
      }else{
      
      dP = I * Aeff(gain,nu/GHz) * StepNu * Calib_Dec2012(nu/GHz);
     
    }
    
    if( antenna == 3) {
      //cout << dP << endl;
      
      //cout <<l << "\t"<< NuMin <<  "\t" << StepNu/GHz << "\t"  <<  nu << "\t gain: " << gain << "\t Aeff: "<< Aeff(gain,nu/GHz) << "\t calib: "<<  Calib_Dec2012(nu/GHz) << "\t" << dP<< endl;


      hh1_exp->SetBinContent(l+1, dP              );
      hh2_exp->SetBinContent(l+1, dP / Calib_Dec2012(nu/GHz) );
      //hh1_exp->Fill(nu/GHz,10.*log10( dP*1e3 )              );
      //hh2_exp->Fill(nu/GHz,10.*log10( dP*1e3 / Calib_Dec2012(nu/GHz) ) );
    }

    P_cal = P_cal + dP;
    P     = P     + dP/Calib_Dec2012(nu/GHz);
  }

  P_cal_dbm = 10.*log10( P_cal*1e3 );
  V_cal = sqrt( P_cal * R );

  P_dbm = 10.*log10( P*1e3 );
  V = sqrt( P * R );

  //cout << V_cal <<  endl;
  
  return (V_cal);

}

double GetCalibIntegral(int selfreq, vector<double> MinFreq, vector<double> MaxFreq){
  
  // frequencies in GHz
  double bandwidth = 10. *1e-3; // GHz
  double NuMin = 2.0 - bandwidth/2.; //GHz
  double NuMax = 20.0 + bandwidth/2.; //GHz
  int NNu = (NuMax - NuMin)/bandwidth;  
  double StepNu = bandwidth * GHz;  //Hz
  
  double P_int = 0;
  for (int l=0;l<NNu;l++) {
    
    double nu = (NuMin + l * StepNu/GHz) * GHz; // nu in Hz
    double gain = tgAntenna->Eval(nu); //nu in Hz
    
    bool cut=false;
    for (unsigned int f=0; f<MinFreq.size();f++){
      if(nu/GHz>MinFreq.at(f) && nu/GHz<MaxFreq.at(f))  //nu in GHz
	cut = true;
    }
    
    double dP=0;
    if(selfreq==1 && cut==true)
      dP=0;
    else
      dP = Aeff(gain,nu/GHz) * StepNu * Calib_Dec2012(nu/GHz);
     
    //cout <<l << "\t"<< NuMin <<  "\t" << StepNu/GHz << "\t"  <<  nu << "\t gain: " << gain << "\t Aeff: "<< Aeff(gain,nu/GHz) << "\t calib: "<<  Calib_Dec2012(nu/GHz) << "\t" << dP<< endl;
    
    P_int = P_int + dP;
  }

  return P_int;

}



